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USERS MANUAL FOR THE NASA LEWIS THREE-DIMENSIONAL 
ICE ACCRETION CODE (LEWICE3D) 


Colin S. Bid well 
Mark G. Potapzcuk 

National Aeronautics and Space Administration 
Lewis Research Center 
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SUMMARY 


A description of the methodology, the algorithms, and the input and output data along with 
an example case, for the NASA Lewis three-dimensional ice accretion code (LEWICE3D) has 
been produced. The manual has been designed to help the user understand the capabilities, the 
methodologies and the use of the code. 

The LEWICE3D code is a conglomeration of several codes for the puipose of calculating 
ice shapes on three-dimensional external surfaces. A three-dimensional external flow panel code 
is incorporated which has the capability of calculating flow about arbitrary three-dimensional lift- 
ing and nonlifting bodies with external flow. A 4th order Runge-Kutta integration scheme is used 
to calculate arbitrary streamlines. An Adams type predictor-corrector trajectory integration 
scheme has been included to calculate arbitrary trajectories. Schemes for calculating tangent tra- 
jectories, collection efficiencies and concentration factors for arbitrary regions of interest for sin- 
gle droplets or droplet distributions have been incorporated. A heat transfer algorithm based on 
the NASA Lewis two-dimensional ice accretion code (LEWICE) can be used to calculate ice 
accretions along surface streamlines. A geometry modification scheme is incorporated which cal- 
culates the new geometry based on the ice accretions generated at each section of interest. 

The three-dimensional ice accretion calculation is based on the two-dimensional LEWICE 
calculation. Both codes calculate the flow, pressure distribution, and collection efficiency distribu- 
tion along surface streamlines. For both codes the heat transfer calculation is divided into two 
regions, one above the stagnation point and one below the stagnation point, and solved for each 
region assuming a flat plate with pressure distribution. Water is assumed to follow the surface 
streamlines, hence starting at the stagnation zone any water that is not frozen out at a control vol- 
ume is assumed to run back into the next control volume. After the amount of frozen water at each 
control volume has been calculated the geometry is modified by adding the ice at each control vol- 
ume in the surface normal direction. 
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SYMBOLS 


A, A2 

A 0 

A SU r 

B, B2 
B k 

Cp 

DCOR 

DFINE 

DS 0 

DS m 

DICE 

DICES 

dice 

f 

h c 

i 

IFLOW 

IRUN 

ICE 

IMOD 

ISUP 

ISLO 

ISTRF 

IST1, IST2 

Lf 

Ly 

m 

NBR 

NBC 

NPSEC 

NPTS1, NPTS2 

PIN, PIN 1, 

P1JP2 

PLN 

q c 

T 

TL 


Parametric slope matrix for a line in space 
Trajectory flux tube area at surface 
Trajectory flux tube area in free stream 
Surface area of a segment on the streamline 
Parametric intercept matrix for a line in space 
Vortex pair strength for kth lifting strip 
Specific heat, J/kg 

Coarse step size for tangent trajectory search 
Fine step size for tangent trajectory search 
Trajectory flux tube width in free stream 
Trajectory flux tube width at surface 
Ice thickness array along streamline 
Ice thickness save array for streamlines 
Ice thickness array 
Freezing fraction at a segment 
Convective heat transfer coefficient 
Enthalpy 

Flow field calculation flag 
Trajectory calculation control flag 
Ice accretion calculation control flag 
Geometry modification control flag 
Upper streamline release point array counter 
Lower streamline release point array counter 
Streamline calculation control flag 

Rags denoting two closest sections of interest to a given N- 
line 

Heat of fusion, J/kg 
Heat of vaporization J/kg 
mass flow rate 

Number of rows of trajectories to be released at a section of 
interest 

Number of columns of trajectories to be released at a section 
of interest 

Rag describing region of interest 

Number of points in two the two closest streamlines to a giv- 
en N-line 
Point arrays 

Array containing coefficients of equation of a plane 
Convective heat flux, W/m 2 
Temperature, K 

Parametric distance of an N-line from two closest sections of 
interest 
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SYMBOLS CONTINUED 


u 


V 

VCRIT 

XNEW, 

YNEW, 

ZNEW, 

XSCI, 

YSCI, 

ZSCI 

XTIP, 

YTIP, 

XTIP 

XSEC, 

YSEC, 

ZSEC 

XSTOP 

XL,YL, ZL 

XN,YN,ZN 

X i 

x ci 

x si 

x j 

Xi 

x s 

a 

P 

Pi 

As 

At 

V 
v 


Subscripts 


a 

aw 

c 

e 


Weighting factor in j direction for collection efficiency inter- 
polation 

Weighting factor in i direction for collection efficiency inter- 
polation 
Velocity M/S 

Stagnation point search velocity criteria 
Coordinate arrays describing the on-body streamline. 


Arrays describing upstream release points 

for trajectories passing through points describing the 

section of interest. 

Arrays describing upstream release points for tangent 
trajectories 

Arrays describing the region of interest. 


Stream wise stopping point for all trajectory and streamline 
calculations 

Scan line arrays for stagnation point search 
Arrays of surface normals for streamline 
Surface coordinates of impingement cell 
Location of centroid of A m 
Location along surface streamline 
Displacement vector from x(ij) to x(i,j+l) 

Displacement vector from x(ij) to x(i+l,j) 

Displacement vector from x(i j) to x s 

Pitch angle of geometry 

Collection efficiency 

Density of ice at segment i 

Length of segment along surface streamline 

Time increment for ice accretion 

Sideslip angle of geometry 

Rotation angle of surface droplet flux tube 

Source strength of a panel 


Air 

Adiabatic wall 
Critical; convection 

Evaporation: condition at the edge of the boundary layer 
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SYMBOLS CONCLUDED 


i Ice 

(i) Control volume 

r in Runback into control volume 

r out Runback out of control volume 

sur surface condition 

s Static condition 

T Total condition 


I. INTRODUCTION 

The LEWICE3D code is a conglomeration of several codes for the purpose of calculating 
ice shapes on two-dimensional external surfaces. A three-dimensional external flow panel code is 
incorporated which has the capability of calculating flow about arbitrary three-dimensional lifting 
and nonlifting bodies with external flow. A 4th order Runge-Kutta integration scheme is used to 
calculate arbitrary streamlines. An Adams type predictor-corrector trajectory integration scheme 
has been included to calculate arbitrary trajectories. Schemes for calculating tangent trajectories, 
collection efficiencies and concentration factors for arbitrary regions of interest for single droplets 
or droplet distributions have been incorporated. A heat transfer algorithm based on the NASA 
Lewis two-dimensional ice accretion code (LEWICE) can be used to calculate ice accretions 
along surface streamlines. A geometry modification scheme is incorporated which calculates the 
new geometry based on the ice accretions generated at each section of interest. 

The three-dimensional ice accretion calculation is based on the two-dimensional LEWICE 
calculation. Both codes calculate the flow, pressure distribution, and collection efficiency distribu- 
tion along surface streamlines. For both codes the heat transfer calculation is divided into two 
regions, one above the stagnation point and one below the stagnation point, and solved for each 
region assuming a flat plate with pressure distribution. Water is assumed to follow the surface 
streamlines, hence starting at the stagnation zone any water that is not frozen out at a control vol- 
ume is assumed to run back into the next control volume. After the amount of frozen water at each 
control volume has been calculated the geometry is modified by adding the ice at each control vol- 
ume in the surface normal direction. 

The basic methodology of the three-dimensional ice accretion analysis is to divide the 
three-dimensional ice accretion process into two-dimensional processes along streamlines of inter- 
est (fig. 1). The user inputs regions of interest on the three-dimensional body (e.g. leading edge 
points). A streamline is then calculated along the body’s surface from the centroid of this region of 
interest. Impingement rates and velocities are calculated along this streamline. This information is 
input to a two-dimensional heat transfer module which calculates ice growth along the streamline. 
This information is used to generate a new geometry at the streamline location. This process is 
repeated for each streamline of interest on the three-dimensional body. Upon completing the ice 
growth calculations the geometry is modified and the flow field is updated. The above steps are 
repeated for as for each time step. 
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(a) Clean swept airfoil with three streamlines (top view). 



(b) Clean swept airfoil with three iced streamlines (top view). 



(c) Iced swept airfoil resulting from three iced streamlines. 
Figure 1. - Airfoil with several sections of interest. 


The three-dimensional (3D) analysis then, can be broken into 6 basic steps. First, a flow 
field is generated for the body. Second, impingement efficiency is calculated at the region of inter- 
est. Third, a streamline is calculated at the region of interest. Fourth, impingement rates along the 
streamline of interest are found by interpolation. During the Fifth step ice accretion along the 
streamline is calculated using the two-dimensional (2D) heat transfer module. The Sixth step 
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involves generating a new body from the ice accretion information. 

A 3D Hess-Smith panel code (ref. 1-5) is used to generate the flow field used in the trajec- 
tory and heat transfer calculations. The code can accommodate lifting and non-lifting geometries 
or combinations thereof such as entire airplanes (fig. 2). If desired, a Prandtl-Glauert correction can 
be made for compressible cases. The code can also handle leaking panels to emulate inlets for 
instrument orifices. The code also has a variable dimension feature which allows easy adaption to 
different computers or problems. 



(a) Finite swept wing 



(b) Twin Otter aircraft 

Figure 2. - Panel representation of different types of geometries 


The trajectory code is basically that developed by Hillyer Norment (ref. 5) with one addi- 
tional feature. The code uses the Hess-Smith flow field along with an Adams-type predictor-cor- 
rector algori thm developed by Krogh (ref. 6). An added feature is the ability to calculate local 
collection efficiency from the impacting trajectories. The code is used here to generate an array of 
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impingement efficiencies for each region of interest. 

The surface streamline is calculated using a 4^* order Runge-Kutta integration scheme. The 
streamline integration is carried forward from the stagnation region for both the upper and lower 
surfaces at the region of interest. 

A linear interpolation scheme is used to determine the collection efficiency along the 
streamline from the matrix of collection efficiencies generated above in the trajectory step. 

The 2D ice accretion calculation is basically that of the LEWICE program generated at 
Lewis. This code is described in detail in reference 7. 

The new geometry is generated from ice accretion information and from the surface normal 
information and final trajectory angle information. Each new point on the streamline is generated 
by adding the ice accretion multiplied by either the surface normal vector or by the final trajectory 
tangent vector to the old streamline point. 

II. PROGRAM UNITS 

A. Introduction 

There are nine basic program units comprising the 3D ice accretion calculation: READIN, 
FLOW, SETFLO, BETAC, STREM3D, STREM2D, BSTREM, LEWICE2D,AND BODMOD. A 
brief description of each of these modules is given along with a flow chart. Figure 3 shows an over- 
view of the LEWICE3D job stream. Section J contains tables giving a brief description of each sub- 
routine used in the above modules. 



Figure 3. - LEWICE3D segmentation tree structure. 













B. Subroutine READIN 


1. General Discusion 

The module READIN reads the job control file (unit INPUT) and initializes important pro- 
gram control variables. All input data (unit INPUT) is in “NAMELIST 5 format. Three “NAMEL- 
ISTS” are input from unit INPUT: IMPING, TRAJ and ICEIN. A brief description of the variables 
in each of the NAMELISTS is given in the INPUT FILE section. IMPING contains control vari- 
ables for the Hillyer Norment Trajectory codes. These variables are described in the user’s manual 
for the trajectory codes (ref. 5). TRAJ contains the control variables for the overall calculation 
including how many stations are to be used, number of trajectories to be used, whether to run the 
flow field code, LEWICE or streamline calculations, etc. ICEIN contains control variables for the 
2D LEWICE calculation. These variables are described in the LEWICE manual. 

2. Printed Output 

Subroutine READIN prints job control information to several output files (OUPUT, JOB- 
SUM). This information is self explanatory. 

C. Subroutine FLOW 

1. General Discussion 

Subroutine FLOW is essentially the HESS-Smith 3D panel code put into subroutine form. 
Hillyer Norment gives a good description of the Hess-Smith code in his user’s manual (ref. 5, pages 
10-14), and this description is repeated here for completeness. Figure 4 shows a flow chart of the 
flow field (subroutine FLOW) and velocity calculations (subroutine FLOVE2 and FLOVEL). Sub- 
routine FLOVE2 is the original velocity calculation alogorithm developed by Hillyer Norment (ref. 
5). Subroutine FLOVEL is a vectorized version of FLOVE2 which was developed at LEWIS by 
Bidwell and Mohler. Subroutine FLOVEL evaluates velocities about 20 % faster than FLOVE2 on 
computers that do not support vectorization and about 80% faster on machines supporting vector- 
ization. The algorithm cannot be used in cases where the piecewise linear vorticity option has been 
chosen (i.e. PESWIS = TRUE). 
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(a) Subroutine FLOW 



(b) Subroutine FLOVE2 

Figure 4. Flow field and velocity calculation segmentation tree structures. 
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(c) Subroutine FLOVEL 

Figure 4. -Concluded. Flow field and velocity calculation segmentation tree structures. 


The methods and codes of Hess (ref. l)and Hess and Smith (ref. 2,3) are used for calcula- 
tion of lifting and nonlifting potential flow about arbitrary three-dimensional bodies. Lifting bodies 
(i.e., airfoils) alone, nonlifting bodies alone, or combinations of lifting bodies with nonlifting bod- 
ies (e.g., combinations or airfoils and fuselages) can be treated. Effects of flow into an inlet, for 
example an instrument aperture, can be accounted for provided the intake flow rate, in terms of 
fraction of free stream air speed, is specified. The method is restricted to subsonic airspeeds, but 
for free stream Mach numbers greater than 0.5, the Prandtl-Glauert method is used to correct 
approximately for compressibility effects. Since potential flow is computed, neither viscous effects 
nor turbulence are treated. 

The code requires input of a digital description of the body surface, and for purposes of 
organizing the data as well as for computing flow, the body surface is partitioned into sections 
which are designated as either lifting or nonlifting. In either case, the surface is represented by con- 
tiguous, plane quadrilateral panels, usually called elements (fig.2). For nonlifting sections there are 
few restrictions on the manner in which the elements can be arranged to represent the surface other 
than those required for organization. Lifting sections are restricted as follows: each must consist of 
strips of elements, the strips being oriented parallel to the chondwise direction of the airfoil each 
strip must have the same number of elements and wake elements must be included after the trailing 
edge of each strip. Both lifting and nonlifting portions of the body may be described by more than 
one section. 

Each on-body element (which is in the flow) is taken to be a potential flow source. The 
source is a distributed one, with the distribution being uniform over the surface of the element, and 
each element, for example the j** 1 , is characterized by a unique source density, Oj. In addition, each 
strip of elements in a lifting section is characterized by having a unique value of lift vorticity asso- 
ciated with it. This quantity, for example for the K* lifting strip, B®, represents vortex strength 
per unit path length around the strip (fig. 5), and it represent the sum of contributions from all pan- 
els in the strip. Velocities induced by these vorticities are treated as onset flows. Thus, there is an 
onset flow from each lifting strip plus the free stream onset flow. It is necessary to compute an inde- 
pendent source density for each of these onset flows for each on-body quadrilateral panel: if there 
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are N on-body panels, K lifting strips, and one free stream flow, N(K +1) values of o must be com- 
puted. Source densities are determined by solving large systems of linear equations that represent 
the effects of all onset flows on all panels, plus the mutual interactions of all distributed sources, 
under boundary conditions for zero flux through the centroid (also called control point) of each on- 
body panel, or specified fraction of free stream flux through each inlet panel. 


n 



Figure 5. - Organization of m and n lines in a lifting section. A lifting strip is delineated by 
sequential n lines, and extends over the complete circuit from m = 1 at the trailing 
edge, along the underside to and around the leading edge, back to the trailing edge, 
and finally back to the furthest aftward extent of the wake. 


Determination of vortex strengths requires an additional constraint, the Kutta condition, 
and this is supplied by user-selection of one of two optional methods which are designated as “flow 
tangency” and “pressure equality.” 

Lift vorticity is computed by a novel method developed by Hess (ref. 1). To circumvent 
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problems that have been found to result from use of vortex filaments in prior work, and to ensure 
that potential flow results from the vorticity distribution and that individual infinitesimal vortex 
lines either form closed curves or go to infinity, Hess has developed a method by which vortex 
sheets on the body and wake surfaces can be expressed in terms of dipole sheets on the same sur- 
faces. Hess summarizes the method as follows: 

“A variable- strength dipole sheet is equivalent to the sum of: (1) a variable- 
strength vortex sheet on the same surface as the dipole sheet whose vorticity has a 
direction at right angles to the gradient of the dipole strength and a magnitude 
equal to the magnitude of this gradient, and (2) a concentrated vortex filament 
around the edge of the sheet whose strength is everywhere equal to the local edge 
value of dipole strength.” 

Mathematical details are given in appendix A of reference 1. 

For particular body geometry and orientation relative to the free stream, the source densities 
and vortex strengths are calculated only once, and then these can be used to calculate flow velocity 
at any space point exterior to the body. The primary functions of the DUGLFT codes are to calcu- 
late the Oj and and store these quantities, along with other requisite data, for use by subroutine 
FLOVEL in calculating flow velocities. Subroutine FLOVEL is called as needed by programs 
TRAJEC, CONFAC, and ARYTRJ to provide flow velocities for trajectory and flow velocity array 
calculations. 

In calculating each flow velocity, contributions from all quadrilateral elements are summed. 
There are three sets of algorithms for computing contributions from individual elements: (1) for 
elements that are close to the calculation point, detailed calculations are used that account for exact 
element geometries, (2) for elements at intermediate distances multipole expansions are used, and 
(3) for remote elements the point source approximation is used. Mathematical details are given in 
references 1,2 and 3 with emphasis on lifting flow in reference 1 and emphasis on nonlifting flow 
in references 2 and 3. The reader is strongly urged to study these references closely before attempt- 
ing to use this code. Reference 4 consists of a code users manual for the lifting flow calculations 
described in reference 1. 

Calculation accuracy is discussed in the Validation section in the Hillyer manual (ref. 5). 
Of course accuracy also depends on the fineness of resolution of the element description of the 
body, and naturally some compromise is called for. The smaller the elements the finer the resolu- 
tion, and the fewer of them for which the most exacting of the three algorithms must be used. On 
the other hand, the number of elements increases inversely as the square of their linear size. In past 
studies on airplanes we have used the following paneling criteria For those parts of the airplane 
traversed by particle trajectories, we try to keep the element edges between 6” and 8” in length. 
Where allowed by simplicity of surface shape, remote elements can be larger. Remote downstream 
complexities of shape are ignored or treated approximately. For example, if interest is confined to 
the forward fuselage, then the remainder of the fiiselage can be represented as a cylinder of constant 
cross-section which is extended to approximately five time the length of the of the nose section (as 
recommended by Hess and Smith, ref. 2), and the wings can be ignored entirely. The following are 
basic requirements of the method that apply to all calculations 
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1. A uniform, unit-speed free stream approximately in the direction of the positive x-axis. 

2. Normalization of all velocities to be consistent with the unit free stream speed. 

3. Normalization of all distances by a user-specified characteristic dimension of the body. 

Surface point coordinates may be recorded in any convenient units and can be appropriately trans- 
lated and scaled, to meet requirement 3 above, during processing via use of SR’s PATPRS and 
DATPRS. These subroutines also allow rotation of the body about the y axis to adjust attitude 
angle. The coordinate system used for the calculations is described on pp. 19-20 (ref. 5). 

The unit free stream speed is assumed by program DUGLFT, and the distance normaliza- 
tion, if required, is done during preliminary data processing as indicated above. For trajectory cal- 
culations, the user specifies the true free stream speed and the normalization length, and the codes 
automatically handle any additional normalizing or scaling that is required. 

The module FLOW computes a flow field for the geometry input on unit NGEOM (DUG- 
LIFT format) and saves it on unit FLOWF. This module is executed only if EFLOW=l (NAMEL- 
IST TRAJ). If the flow field code is not executed (i.e. IFLQW=0) then the flow field must be 
provided on unit FLOWF. The Hess-Smith 3D flow field code is used in subroutine form here to 
generate the flow field. The execution time for the flow field calculation is proportional to the 
square of the number of panels. A 1200 lifting panel model with one section required 80 seconds 
of CPU time on the CRAY XMP while that for a 3200 panel model with one section required 630 
seconds. Two basic requirements have been found to date for the ice accretion calculations. The 
first is a numerical one for the panel method used. The requirement is that the aspect ratio of any 
panel should not be greater than 100. The source strength calculation will converge for larger 
aspect ratios but the vortex calculation will not. The requirement grows more stringent with angle 
of attack and ratios as high as 100 may not be allowed for some geometries. The second require- 
ments is that to produce a smooth beta curve there must be approximately one panel per trajectory 
released in the z-direction (i.e. if 20 trajectories were to be released between the impingement lim- 
its then 20 panels are required between the impingement limits at the surface to ensure a smooth 
beta curve. 

2. Printed Output 

Subroutine flow produces several output files (units FLSUM, FLOWF). Unit FLSUM is a 
summary of the flow field computation and contains varying amounts of information depending on 
the flags set in the flow input file (unit NGEOM). Any error messages from the flow field compu- 
tation will be found in unit FLSUM. Unit FLOWF contains the flow field information in binary 
format to be used in the calculation of velocities in the trajectory routines. A description of these 
files is taken from reference 5. 

The flow field calculation summary output (unit FLSUM) consists of two main parts, plus 
a s ummar y of input control data, various error condition messages, and optional outputs of data 
that are used for debugging. 
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1. The first printout is a summary of input control data, and is self-explanatory. 

2. Next, which is the first main part, is a printout (from NOLEFT and LIFT) of element 
data. Elements are designated as lifting or non-lifting. 

A short table follows (from INPUTL) titled TABLE OF INPUT INFORMATION, 
which s ummari zes the data in terms of section type, number of elements per section, 
number of strips, etc. 

3. In the course of computing velocities induced by each element on all others, additional 
su mmar y information is printed (from VDMX) for each section. For lifting strips, this 
includes information on ignored elements which does not appear elsewhere. 

4. If the piecewise linear option for determination of spanwise variation of vortex strength 
is used, strip widths, W^, and parameters D^, E^, (ref. 3, sec. 7.11) are printed for 
each strip for each section, along with a summary of edge conditions (NLENE1 and 
NLINEN). 

5 A short statement is printed (from COLSOL) regarding the dimensions of the matrix 
that are solved to determine element source strengths, and the number of right-hand- 
sides (i.e., number of uniform onset flows plus number of lifting strips) for which the 
solutions are obtained. 

6. The second main printout (from PRINTL) contains the final results of the calculations. 
A printout for each on-body element is labeled as follows. 


XO, YO, ZO 

Control point coordinates. 

VX, VY, VZ 

Flow velocity components at the control point 

VT 

Velocity magnitude 

VTSQ 

Square of velocity magnitude 

CP 

Pressure coefficient = 1.0 - VTSQ 

DCX, DCY, DCZ 

Direction cosines of the velocity components 

NX, NY, NZ 

Components of the unit normal to the plane of the el- 
ement 

SIG 

Source density 

VN 

Velocity component in the direction of the unit nor- 
mal 
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AREA 


Area of the element 


Printouts for off-body and Kutta points are similarly labeled. 

Also printed are vector components for pressure force and moment for each strip, each 
section and for the entire body, as well as a table of vortex strength per unit length, B , 
for each lifting strip. 

7. Error messages (ref. 4). 

(a) Message: MISMATCH OF ELEMENTS IN A LIFTING 

STRIP IS DETECTED. ELEMENTS FORMED 
XXX, ELEMENTS INPUT XXX, COMPUTATION 
TERMINATED. (SR INPUTL) 

Cause of error: Inconsistent input data. The program sums the num- 

ber of on-body elements plus the wake elements 
specified on card 8. This sum does not match with the 
elements formed from the input coordinates. 

Action: Check the lifting body information card (card 8) and 

the quadrilateral comer point coordinates cards 
(cards 12). The number of points on an n-line should 
equal the number of elements plus 1. 

For example: If in a lifting section each lifting strip 
consist of 10 on-body elements and 1 wake element, 
the total number of elements is 11, and there should 
be 12 points on each n-line input via cards no. 12. 

(b) Message: ERROR IN IGNORED ELEMENT COUNT XXX, 

SHOULD BE XXX. (SR LIFT) 

Cause of error: Erroneous specification of ignored element informa- 

tion. 

Action: Check card 10 to make sure the ignored element in- 

formation is properly specified. 

(c) Message: LABEL ERROR IN NONLIFTING VFORM. (SR 

VFMNLF) 

LABEL ERROR IN LIFTING VFORM (SR VFM- 
LFT) 

Cause of error: Geometric data for each element strip, preceded by a 
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lifting or nonlifting label are stored on unit 4. The er- 
ror occurs when a labeling mix-up is detected during 
input of the data from unit 4 for calculation of veloc- 
ities. That is, data for a strip labeled lifting are en- 
countered during computation for a nonlifting 
section, or vice versa. 

Action: Check that the number of lifting strips specified on 

card no. 8 for each lifting section corresponds with 
the cards no. 12 input. 

(d) The following messages pertain to errors in specification of variable 
dimensions (SR CKARRY). 

ELEMENT CAPACITY, NONX = XXX IS LESS THAN TOTAL 
ELEMENTS, NON= XXX 

STRIP CAPACITY, NSTX = XXX IS LESS THAN TOTAL STRIPS, 
NSTRP = XXX 


LIFTING SECTION CAPACITY, LFSX= XXX IS LESS THAN TOTAL 
SECTIONS, ISECT = XXX 

LIFTING STRIP CAPACITY, NOBX = XXX IS LESS THAN TOTAL 
LIFTING STRIPS, LSTRP = XXX 

N2BX = XXX IS NOT GE TWICE NOBX AS REQUIRED, NOBX 

= xxx 

NSLX = XXX IS LESS THATN THE MAX. NO. OF STRIPS IN A 
LIFTING SECTION, WHICH IS XXX 

CAPACITY OF ARRAY WKAREA, NWAX = XXX, USED BY COLSOL 
TO DETERMINE SOURCE STRENGTHS, IS INSUFFICIENT. IT 
MUST BE GREATER OR EQUAL TO XXX 

NWAX = XXX IS NOT GREATER OR EQUAL TO NO. OF LIFTING 
STRIPS = XXX CUBED, AS REQURED FOR THE PRESSURE 
EQUALITY KUTTA OPTION. 

Cause of error: Array dimensions are inadequate to accommodate 

the input data. 

Action: Check array dimensions and variable array parame- 

ters against the storage demands of the element data 
input via cards no. 12. Also check input parameter 
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(e) Messages: 


Cause of error: 


Action: 


(f) Message: 


Cause of error: 


Action: 


(g) Message: 


Cause of error: 


Action: 
(h) Message: 


LIFSEC, NSORCE, NWAKE, NSTRIP, and IX- 
FLAG. 

XXX ANGLES OF ATTACK HAVE BEEN SPECI- 
FIED, ONLY ONE IS ALLOWED SINCE COM- 
PRESSION EFFECTS ARE CONSIDERED. 

ANLGE OF ATTACK, + - XXX, + - XXX, +- XXX- 
IS INAPPROPRIATE FOR A CASE WITH COM- 
PRESSION CORRECTION. (SR CKARRY) 

Only one uniform onset flow (i.e. free stream) is al- 
lowed if the compressibility correction is applied 
(MACH > 0.0 on card no. 2). Moreover, the direction 
cosines (ALPHAX, ALPHAY, ALPHAZ) of this on- 
set flow must be (1.0, 0.0, 0.0; card no. 4). 

Set IATACK = 1 on card 2, and/or specify direction 
cosines on card 4 as stated above. 

THE NUMBER OF KUTTA POINTS SPECIFIED 
IS INCORRECT AND SHOULD BE XXX (SR CK- 
ARRY) 

The flow tangency Kutta option has been specified, 
and the number of Kutta points specified by input 
(cards no. 9, 13 and 14) does not equal the number of 
lifting strips. 

Check parameter KUTTA on card 9, and the number 
of KUTTA data points on cards 1 3 and 1 4, against the 
number of lifting strips input via cards no. 12 (Do not 
count extra strips.) 

ERROR IN VFORM. THE ELEMENTS FORMED 
DO NOT CORRESPOND TO THE NO. OF BODY 
ELEMENTS. (SRS VFMNLF AND VFMLFT) 

Element tally recorded by SR INPTL does not match 
with tally recorded from input of data from unit 4 
during velocity calculation. 

Check lifting strip specification data on card 8 for 
consistency with cards no. 12 input data. 

AFTER XXX INTERATIONS, DELTA B STILL 
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Cause of error: 

Action: 

(i) Message: 


Cause of error: 

Action: 

(j) Message: 

Cause of error: 

Action: 

(k) Message: 

Cause of error: 

Action: 

(l) Message: 


DID NOT CONVERGE TO THE GIVEN CRITERI- 
ON/ LARGEST DELTA B— +- XXX .XXX/ PRO- 
GRAM PROCEEDS WITH THE MODS 
CURRENT VORTEX STRENGTH. (SR PKUTTA) 


Non-convergence of vortex strengths, B, calculation 
via the pressure equality Kutta condition method 
(ref.4, sec. 7.13.2). 

Check the cards no. 12 input data. 

THIS CODE SHOULD BE APPLIED TO FIRST 
STRIP 


or, 


THIS CODE SHOULD BE APPLIED TO LAST 
STRIP. (SR DKEKFK OR PSONST) 

Improper specification of NLENE1 or NLINEN for 
piecewise linear option. Specifically, either NLINE1 
= 2 or NLINEN =3 is specified, both of which are for- 
bidden. 

Check card 8 specifications 

XXX ON-BODY POINTS MISSED. EXECUTION 
TERMINATED. (SR PRINTL) 

The number of on-body source elements tallied dur- 
ing final printout does not agree with the count tallied 
during input. 

Check input data. 

XXX KUTTA POINTS MISSED, EXECUTION 
TERMINATED. (SR PRINTL) 

The number of Kutta points tallied dining the final 
print out does not agree with the number specified by 
parameter KUTTA on card 9. 

Check the number of Kutta points input via cards 1 3 
and 14 against parameter KUTTA. 

XXX OFF-BODY POINTS MISSED, EXECU- 
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HON TERMINATED. (SR PRINTL) 

Cause of error: The number of ofF-body points tallied during final 

printout does not agree with the number tallied dur- 
ing input 

Action: Check input data. 

8. Optional Printouts for Use in Debugging 

(a) Geometrical data for each element. (IOUT = TRUE, card 3, SR ENPUTL). 
For each nonlifting element is printed the element sequence number and 
twenty-nine geometric quantities (ref.4, sec. 9.51) and for each lifting ele- 
ment is printed the element sequence number and forty-five geometric quan- 
tities (ref. 3, sec. 7.2). 

(b) Source induced velocity matrix, V t j. (MPR = 1, card 2; SR PNTVU) 

COLUMN Matrix column number (j) 

CNTRL PT Control point number (i) 

VXS, VYS, VZS Velocity components 

_»F 

if IjJFSEC is greater than 0 (card 2), dipole induced velocity matrices, V Vj , 
V ik , also are printed. 

STRIP Lifting strip number 

CNTRL PT Control point number 

VXF, VYF, VZF First velocity components 

VXS , VYS , VZS Second velocity components 

(^) (°°) 

(c) Onset flow matrices, Vj , V,- . (MPR = 3, SR UNIFLO) 

ONSET FLOW NO. 

CONTROL POINTS Control point number 

X-FLOW, Y-FLOW Onset flow velocity components 
Z-FLOW 

(d) Dot product matrices, Ay, N^, N^ 00 -* (MPR > = 2, card 2 SR AUMX and 
NIKMX) 
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COLUMN 
A U 

FLOW NO. 
RIJ 


Matrix column number (j) - 
Elements of Ay 
Onset flow number (k) 

Right side of equation (7.12.5) 


(e) Source density matrix (MPR > 2, card 2; PGM SIGMAL) SOLUTION 
OBTAINED AFTER COLSOL FLOW NO. Onset flow number. Element 
source densities, Oj®, o/°°\ are printed eight to a line. 

The following data are stored on unit FLOWF in binary format for use later in the velocity 
calculations. Actual record structures are more easily determined by examining the SR SETFLO 
FORTRAN listing. 


CASE Body identifier (input card 2) 

ISECT Number of sections (lifting plus nonlifting) 


LIFSEC Number of lifting sections (input card 2) 


ALPHAX(l) Uniform onset flow direction cosines (card 4) 
ALPHAY(l) 

ALPHAZ(l) 


SYM1 

SYM2 

NSYM 

NSTRP 

BETAM 

BETSQ 


Floating point equivalent of input parameters NSYM1 and NSYM2 
(card 2) 

Total number of symmetry planes 

Total number of strips, including extra lifting strips if input. 
SQRT(1-N m 2 ) where N^j is free stream Mach number 
1-N m 2 


NLT Number of elements on each strip, including extra strips, and 

(NSTRIP) ignored and wake elements are counted. It is negative for the last 

strip of each section. 


NTYPE Section type indicator. 
(ISECT) 

0 for nonlifting 
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1 for lifting 

NLINE Number of strips in a section, not including extra strips 
(ISECT) 

If LIFSEC GT 0: 

IGW If true, there are ignored elements. 

LASWAK If true, the semi-infinite final wake element option is exercised 

PESWIS If true, the piecewise linear method for computing spanwise varia- 
tion of lift vorticity is used. 

NSTRIP See input card 8 
(LIFSEC) 

NLINE 1 (LIFSEC) 

NLINEN (LIFSEC) 

NSORCE(LEFSEC) 

IXFLAG(LIFSEC) 

IG 1 (I, J) Only if IGW = TRUE (see input card 10) 

IGN(I,J) 

For each nonlifting element, the twenty-nine geometric quantities written on unit 4 
by SR NOLIFT. 


For each lifting element, the fifty-seven geometric quantities written on unit 4 
by SR LIFT. 


Only if the piecewise linear method is used for calculation of spanwise variation of 
vorticity. For each of K strips in J = LIFSEC lifting sections: 

K,(D(U),E(U)JF(IJ),I=1,K) 

where D, E, F are L^, E k , F k of equation 7.11.5 of reference 3. 

KFLOW Number of lifting strips 

KONTRL Number of on-body source elements (not including ignored, wake, 
and extra strip elements) 

COMSIG Combined source densities (ref. 3 eq. 7.13.1) 

(KONTRL) 
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B(KFLOW) Vortex strength per unit length 


3. Peripheral Storage 

In addition to the flow field summary file (unit FLSUM) and the flow field file (unit 
FLOWF) several internal files are needed for the flow field calculation. Subroutine FLOW used 
eleven units for scratch storage. All data stored on these units are in binary format. In the following, 
use of each unit is considered only in terms of the maximum number of data words (numbers) and 
record lengths that would be stored on it. The following variables are defined to aid in this: 

KONTRL Number of quadrilateral elements, not including 

those generated by symmetry, ignored, in the wake 
and in extra strips. 

KUTTA Points defined by input cards no. 13 and 14 at which 

the Kutta condition is to be applied. (KUTTA > 0 
only if the flow tangency option is exercised.) 

NOFF Number of off-body points at which velocity is to be 

calculated as defined by input cards no. 15. 

NON KONTRL + KUTTA + NOFF 

IATACK Number of lifting strips, not counting extra strips, nor 

those generated by symmetry. 

NFLOW KFLOW + IATACK 

Unit 3: NFLOW records each consisting of 3 x NON numbers 

Unit4: There is a record of 29 numbers for each nonlifting quadrilateral element 

plus 

There is a record of 57 numbers for each lifting quadrilateral element (including 

ignored, wake, and extra strip elements) 

plus 

A one word record for each section of elements 
Unit 8 : The larger of 

Two records each of length 3 x NON numbers or 
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NFLOW records each of length 6 x KFLOW numbers or 
KONTRL records each of length KUTTA numbers 


Unit 9: KONTRL records of maximum length KONTRL + 1 numbers 

Unit 10: KONTRL records of maximum length KONTRL + 1 numbers 

Unit 11: 1/2 KONTRL records each of length 3 x NON numbers 

Unit 12: 1/2 KONTRL records each of length 3 x NON numbers 

Unit 13: The larger of 

2 x KONTRL records each of length 3 x NON numbers or 
KONTRL records of maximum length KONTRL + NFLOW +1 numbers 
Unit 14: The larger of 

KFLOW records each of length 3 x NON numbers or 
KONTRL records of maximum length KONTRL +1 numbers 
Unit 1 5 : The larger of 

KFLOW records each of length 3 x NON numbers or 

KONTRL records of maximum length KONTRL + NFLOW +1 numbers 


4. Variable Array Dimensioning 

Subroutine FLOW incorporates variable dimensioning so that the program can be resized 
to fit different sized problems and computers. The calculation and storage of the flow field account 
for almost all storage required by LEWICE3D and hence are the only areas where variable dimen- 
sioning is used. To resize the problem the variables affected by the following dimension sizes must 
be resized in the main program along with the dimension sizes located in the data statement in the 
main program.The following description of the variable array sizes has been taken from reference 
4. Minimum values for the variable dimension parameters are given where these numbers are not 
effected by symmetry. 

LFSX Number of lifting sections 

NL2X NSLX + 2 
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NOBX Total number of lifting strips, not counting extra strips 

NONX Number of on-body elements in the flow (not counting ignored, wake, and 

extra strip elements) plus Kutta points defined by input (flow tangency 
option only) plus off body points (cards 15) 

NSEX Total number of sections (lifting plus nonlifting) 

NSLX Maximum number of lifting strips in any lifting section (including extra 

strips if input) 

NSTX Total number of strips (i.e., n-lines; lifting plus nonlifting plus extra strips) 

NWAX Similar to 2 x (number of on-body quadrilateral elements in the flow (see 
NONX) plus the number of onset flows) 

and 

Cube of the number of lifting strips (not counting extra strips: i.e. 
NOBX**3) if the pressure equality Kutta condition options is selected. 

N2BX 2 x NOBX 

D. Subroutine SETFLO 
1. General Discussion 

The module SETFLO reads the flow field generated by the Hess-Smith code (unit 
FLOWF). This subroutine is always executed and hence a flow field is required on unit (FLOWF). 

1. Printed Output 

SETLFLO output is limited to summary information about the flow field being read. This 
information is self explanatory and is written to units OUTPUT and JOBSUM. 

E. Subroutine BETAC 
1. General Discussion 

The module BETAC drives the trajectory work used in calculating the local collection effi- 
ciency at each station of interest This subroutine is optional (IRUN=2- 10, NAMELIST TRAJ) and 
controls all the trajectory work (fig. 6). This subroutine calculates tangent and impact trajectories 
at the station of interest. The impact trajectory information is then used in the collection efficiency 
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calculation. If module BSTREM or LEWICE2D are to be executed then the collection efficiency 
information is required and BETAC must be executed. 



Figure 6. - BETAC segmentation tree structure. 


a. Search For Upstream Particle Release Points 

BETAC’s first step toward the calculation of collection efficiency is to determine the 
upstream release points for trajectories that pass through the area of interest (fig. 7). There will be 
one release point calculated for each point specified at the section of interest (i.e. if NPSEC=2 then 
2 upstream release points will be calculated; if NPSEC=4, then 4 upstream points will be calcu- 
lated). These upstream release points (XSCI.YSCI, ZSCI) correspond to trajectories that will pass 
through the points defining the section of interest (XSEC, YSEC, ZSEC). The program CONFAC 
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developed by Norment (ref. 5) and put into subroutine form here is used to determine each of the 
upstream points. 
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^ Array defining trajectory release points that pass through points defining region of 

interest: XSC3<ISC t NPSEC(ISC)X YSCI(ISC),NPSEC(ISC)X 
ZSCI(ISC) f NPSEC(ISQ) 

O Array defining release points for tangent trajectories in region of interest: 
XnP(ISC,NPSEC(ISC)) # YTn>aSC,NPSEC(ISC)X ZTIP(ISC t NPSEC(ISC)) 
Array defining impact points for tangent trajectories in region of interest: 
a XTIF(ISC<NPSEC(ISQ) t YTO(KQNPSEC(ISC)), ZTIF(ISC,NPSEC(ISC)) 


Figure 7. - Illustration of tangent trajectory search arrays. 


Program CONFAC is described in detail in reference 5 and will be covered only briefly 
here. CONFAC is an iterative procedure which finds a trajectory that passes within RW*TOL of 
the specified target point (XSEC, YSEC, ZSEC). Subroutine MAP controls the iteration and uses 
subroutine TRAJEC for the calculation of individual trajectories. If convergence has not been 
reached in 50 iterations then CONFAC will continue with the next upstream point. If one or more 
failures occur during the search for each of the upstream points then the program will terminate. 

b. Search For Tangent Trajectories 
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The second step toward calculating the collection efficiency at a section of interest is to 
determine the tangent trajectories. These are limiting trajectories that impact. Trajectories released 
between corresponding upper and lower tangent trajectories will impact the body. Those released 
outside the tangent trajectories will miss the body. There will be one tangent trajectory for every 
point of interest on the body (i.e. if NPSEC=2 then there will be 2 tangent trajectories, an upper 
and a lower. If NPSEC =4, then there will be 4 tangent trajectories, two upper and two lower). The 
tangent trajectories are found by searching for the proper release points along the lines formed by 
the upstream release points XSCI, YSXI, ZSCI. For NPSEC=2 there will be one line. For 
NPSEC=4 there will be 2 lines. Each line is searched in both directions to determine the upper and 
lower tangent trajectories. The tangent trajectory search is handled by the subroutine IMPLIM. 

Subroutine IMPLIM determines tangent trajectories at the section of interest. Subroutine 
IMPLIM is based on the 2D tangent trajectory search routine used in LEWICE (ref. 7). Subroutine 
IMPLIM requires the input of specified lines (in this case the lines are formed by alternating values 
of XSCI YSCI, ZSCI), an initial start point on the specified line (in this case alternating points 
XSCI, YSCI, ZSCI), and the search tolerance DFINE. The algorithm initiates trajectories along the 
specified line midway between the most current “hit” trajectory and the most current miss trajec- 
tory. If the trajectory impacts the body it becomes the current “hit” trajectory. If it misses the body 
it becomes the current “miss” trajectory. The algorithm terminates the search when the distance 
between the initial start points of the current “hit” and “miss” trajectories is less then DFINE. This 
current impacting trajectory is then the tangent trajectory and is denoted by its upstream release 
points XTIP, YTIP, ZTIP. If the subroutine fails to find an impacting trajectory after three steps then 
it will continue on with the next tangent search. If one or more failures occur during the tangent 
trajectory search at the section of interest then the program will terminate. 

c. Calculation Of Impact Trajectories 

The third step is to determine the matrix of release points for the trajectories to be used in 
the collection efficiency calculation. If NPSEC=2, then a string of NBR equi-spaced release points 
will be generated between the upper and lower tangent trajectory release points (XTIP, YTIP, 
ZTIP). If NPSEC=4, then a matrix of NBR x NBC trajectories release points will be determined. 
The rectangle formed by the four upstream tangent trajectory release points (XTIP, YTIP, ZTIP) is 
divided into NBC equi-spaced columns and NBR equi-spaced rows (fig. 8). 
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NBR(ISC) 



YBU(ISC t NPSEC(ISC)X ZBU(ISC,NPSEC(ISC)) 


Array of centroids of impact areas on surface: XBETA(NBR(ISC) - 1 ,NBC(ISC)- 
I), YBETA(NBR(ISC)-l^fBC(ISC)-lX ZBETACNTBROSQ-l^NBCCISC)-!) 


where: 


NBR(ISC) is the number of rows of trajectories to be released in 
impingement array for the ISC* section, NBC(ISC) is the number of 
columns of trajectories to be released in impingement array for the 
ISC* section. 

Figure 8. - Illustration of starting point and impaction point arrays. 


The forth step involves calculating the trajectories for each of the release points generated 
in the above step. Subroutine ARYTRJ (ref. 5) is used to calculate the individual trajectories. The 
impact points corresponding to the release points are stored in arrays XBETA, YBETA, ZBETA 
for use in the calculation of collection efficiency. 

d. Calculation Of Collection Efficiency 

The fourth and final step involves the calculation of collection efficiency at the section of 
interest. Subroutine BETAC calculates the local collection efficiency in two different ways depend- 
ing on the variable NPSEC. The first method (NPSEC = 4) uses a full 3D calculation for which a 
matri x of impact trajectories (NBC x NBR) is required. This method is to be used for areas where 
fully 3D flow is expected. The second method (NPSEC = 2) uses a 2D method in which a single 
string of NBR impact trajectories is required. This 2D method saves considerable computational 
time and is justifiable for cases where small spanwise variations in the flow field are expected 
throughout the section of interest. 
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The full 3D collection efficiency calculations are straightforward and are analogous to 
those of the 2D problem. 3D collection efficiency is defined as the ratio of the particle flux at the 
surface to the particle flux in the free stream. Or if we follow a group of particle trajectories to the 
surface, then the 3D collection efficiency is the ratio of the surface area to free stream area mapped 
out by the particles. 


0<*c> s V^. 


Eq. 1 


If the flux tube consists of four adjacent trajectories in the release matrix (fig.9) then the collection 
efficiency at the surface can be written where corrections for angle-of-attack and yaw have been 
made. 


P (ij) = (cosT' • cosa • A 0 ( ij ) ) / ( A m ( i,j ) ) Eq. 2 


where the collection efficiency is said to be located at the center of the impact region mapped out 
by the four impacting particles. The angles 'P and a refer to the sideslip and pitch angle of the 
geometry relative to the free stream flow vector (fig. 9).The areas A(i j) and A m (i,j) are calculated 
using subroutine AREAP which calculates the area of an arbitrary polygon by dividing it into tri- 
angles and summing the areas of the individual triangles. 
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(a) Side view. 



(b) Top view. 

Figure 9. - Pitch and sideslip corrections for 3D collection efficiency calculation. 


Calculation of collection efficiency using the “pseudo 2D” methods is similar to the 2D 
methods with several corrections. Corrections for angle-of-attack, yaw angle, sweep angle of sur- 
face and deformation of the flux tube are necessary (fig. 10). The resulting collection efficiency 
equation is then 

|3(0 = (cosy- cosa- cosv ■ DS 0 (i))/ (DS m (i)) Eq. 3 

where the collection efficiency is said to occur at the center of the impacting points of the two tra- 
jectories. The distances DS 0 (i) and DS m (i) represent the distance between the upstream release 
point and the impact points respectively for two adjacent trajectories (fig. 10). As for the 3D case 
the angles a and \P represent the pitch and sideslip angles, respectively. Angle v represents the rota- 
tion angle of the droplet trajectory pair relative to the sweep of the leading edge of the airfoil. 
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(c) Top view 


Figure 10. - Corrections for pitch, sideslip, and rotation for “quasi 2D” collection ef- 
ficiency calculation. 


2. Printed Output 

BETAC output is written to several files and includes summary information about the tra- 
jectory calculation, trajectory coordinates, and several error messages. 

Trajectory s ummar y information includes various information about the CONFAC, 
TANTRA, and ARYTRJ runs and is written to units JOBSUM and OUTPUT. The information 
written to the files is self explanatory. 
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Trajectory coordinates along with other pertinent trajectory parameters are written in binary 
format to unit TEMP25. The information is output according the output time step TPRINT and is 
in the following format (SR’s CONFAC, TANTRA, and ARYTRJ) 

IREC,(T (I) J* 1 (I) ,P3(I) ,P5 (I),P2(I) ,P4(I) ,P6(I), VX, V Y, VZ,H,R, AC) I=1,IREC) 


where: 


IREC 

m 

P1(I),P3(I),P5(I) 

P2(I),P4(I)T > 6(I) 

VX,VY,VZ 

H 

R 

AC 


Total number of cards output 

Integration time at the Ith output step 

X,Y, Z components of the particle velocity respectively 

X,Y,Z components of the particle acceleration respectively 

X,Y,Z components of the flow field velocity respectively 

Integration step size at the Ith output step 

Reynolds number at the Ith output step 

Acceleration modulus at the Ith output step 


Several error messages are written to units OUTPUT, JOBSUM. These messages along 
with an explanation and a possible solution are as follows. 

(a) Message: OUTPUT TRAJECTORIES FROM SUBROUTINE 

CONFAC DOES NOT MATCH THE NUMBER 
REQURIED. (SR BETAC) 

Cause of error: Failure to find a trajectory that passes within a given 

tolerance (TOL) for one of the target points at the 
section of interest (XSEC,YSEC, ZSEC). 

Action: Increase error tolerance (TOL), move section of in- 

terests points farther from body (i.e. XSEC, YSEC, 
ZSEC) or increase the trajectory count limit for the 
CONFAC search (ILIM in data statement in SR 
MAP) 

(b) Message: OUTPUT TRAJECTORIES FROM SUBROUTINE 

TANTRA DOES NOT MATCH THE NUMBER 
REQURIED. (SR BETAC) 


Cause of error: Failure to find a tangent trajectory for one of the im- 

pingement limits at the section of interest (X l lP, 
YTIP, ZTIP). 

Action: Increase the tangent trajectory search step sizes 

(DCORJDFINE) or increase die trajectory count lim- 
it for the TANTRA search (KTLIM in data statement 
in SR TANTRA) 
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F. Subroutine STREM3D 


1. General Discussion 

The module STREM3D determines the streamline at the station of interest. The module 
uses a RUNGE-KUTTA integration scheme to calculate the 3D streamlines. Figure 11 shows a 
schematic of the job stream for STREM3D. STREM3D also generates pressure coefficient and 
velocity information at each of the streamlines points. This module is optional (IRUN=1 ,5,7,9: and 
ISTRF=0; NAMELIST TRAJ). If module BSTREM or LEWICE2D ate to be executed then 
streamline information is required and either STREM3D or STREM2D must be executed. 



Figure 11. - STREM3D segmentation tree structure. 


Four steps are involved in determining the 3D streamline: determination of the local stag- 
nation point, integration of the upper surface streamline, integration of the lower surface stream- 
line, and projection of the upper and lower surface streamlines onto the body. 

b. Search For Stagnation Zone 

A marching procedure is used to determine the stagnation zone. This algorithm steps 
towards the body with a step size of H. At each point a vertical scan of velocities is made. A dot 
product is made from consecutive velocities along the scan line. When the dot product reaches a 
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minimum along the scan line, the value is stored and compared to a criterion, VCRIT (currently 
VCRIT =.5). If the value is less than the criterion, then the stagnation point has been found and the 
current scan line points (XL, YL, ZL) and the points where the dot product reached a minimum 
value are stored (ISUP, ISLO). The points ISUP and ISLO are the points in the scan line arrays 
(XL, YL, ZL) where the upper and lower streamlines’ integration will initiate, respectively. If the 
criterion is not met at the current scan line then the procedure steps towards the body and repeats 
the above process. If the algorithm marches to the leading edge without meeting the criterion, 
VCRIT, then a warning message is printed and the stagnation points (ISUP, ISLO) are set to the 
values where the minimum dot product occurs for the current scan line. The above procedure 
essentially searches for the point near the leading edge where the velocity vector diveigence (mea- 
sured by the dot product of consecutive velocities along the scan line) is the greatest. This point 
should be the stagnation zone and streamlines initiated above this point (ISUP) should follow the 
upper surface, while streamlines initiated below this point should follow the lower surface. 

b. Calculation Of Upper and Lower Streamlines 

Once the stagnation point has been found, the upper and lower surface streamlines can be 
determined. Off-body streamlines are used because velocity gradients on the panel edges make the 
integration of on-body streamlines difficult. The streamlines are integrated using a Runge-Kutta 
4th order integration scheme with variable step size from points located on the scan line stored in 
the previous step (i.e. XL, YL, ZL). The upper streamline is integrated from the point ISUP in the 
scan line array. If problems occur at any time during the integration (i.e the streamline penetrates 
the body or the step size is too small) the streamline calculation will be restarted from the next point 
in the scan line array (i.e. ISUP = ISUP +1). The above process is repeated until a streamline is 
integrated to X = XSTOP with no problems. If failure to integrate a “good” streamline occurs for 
ISUP = NPTS (where NPTS is the number of points in the scan line array) then the program will 
terminate. After the upper surface streamline has been found, the lower streamline search begins. 
The lower streamline integration is initiated from the point ISLO in the scan line arrays. As for the 
upper streamline, the lower streamline is integrated until either a X=XSTOP is reached or until a 
problem arises. If a problem occurs then the streamline is restarted from the next point on the scan 
line (i.e. ISLO = ISLO - 1 ). If ISLO drops below 1 during a restart then the procedure will terminate. 

c. Projection Of Streamlines Onto Surface 

The last step is to project the streamlines onto the body. That is, the off-body integrated 
streamlines are projected onto the body in a surface normal direction to produce a streamline which 
has points lying on the panel edges. Projecting the streamline onto the body is done to allow easier 
geometry modifications in later calculations 

To project the streamline onto the body we must first find the portion of the body where the 
streamline is located. This is done by searching for the panel which is closest to the centroid of the 
current section of interest (XC, YC, ZC). The lifting strip that contains the panel is then used for 
the on-body projection. 

Once the local lifting strip has been found, we can construct the on-body points. There will 
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be NSTREM points in the on-body streamline which correspond to one more than the number of 
panels in the lifting strip. Each point on the on-body streamline will lie upon a panel edge connect- 
ing n-lines in the strip. Hence if there are NPTS panels there will be NPTS+1 points in the N-lines 
and hence NPTS + 1 points in the on-body streamline. 

To determine the on-body points corresponding to each panel edge in the strip, we first must 
calculate several parameters for each edge line in the strip: the surface normal of the panel edge 
(XN, YN, ZN), a line containing the panel edge (A, B), and a plane with the normal of the panel 
edge which contains the panel edge (PLN). For the first and last edge line, the surface normal is the 
surface normal for the first and last panel respectively. For internal edges, the surface normal is 
taken as the average of the two surface normals from the panels that form the edge. The line con- 
taining the panel edge is formed from the corresponding points on the panel edge. The edge normal 
plane is determined from the panel edge surface normal and from a panel comer point on the panel 
edge line. 

We can now contsruct the on-body points. Each point on the upper and lower integrated 
streamlines is projected onto the panel edge plane (PLN) using subroutine TRNSF (these points 
are stored in arrays XNEW, YNEW, ZNEW). The algorithm then searches for the two closest 
points in XNEW, YNEW, ZNEW to the line formed by the panel edge (A, B). A line is then formed 
from the two points (A2, B2) and the intersection between the two lines (A, B and A2, B2) is found. 
This point of intersection (PIN) is then the on-body point from the current edge line. This proce- 
dure is repeated for each of the panel edge lines in the lifting strip. 

2. Printed Output 

Output from STREM3D is written to several output files (JOBSUM, OUTPUT) and in- 
cludes summary information about the streamline integration and several error messages. The sum- 
mary information includes coordinates, surface distances (measured from stagnation zone where 
positive surface distance denotes lower surface and negative surface distance denotes upper sur- 
face), pressure coefficients, and surface normals for each point on the streamline. The following 
error messages which can occur are explained with possible solutions. 

Message: STOP IN SUBROUTINE STREAM, ISUP = XXX 

or 

STOP IN SUBROUTINE STREAM, ISLO = XXX 

Cause of error: Failure to find an upper streamline or lower stream- 

line along the line of points (XL, YL, ZL). During the 
search iteration the release point for the streamline 
(ISLO or ISUP) has reached the upper or lower 
bound of the arrays XL, YL, ZL (i.e ISUP = 25 or 
ISLO =1). 

Action: Increase the distance between the body and the sec- 

tion of interest (XSEC, YSEC, ZSEC) 


35 



G. Subroutine STREM2D 


1. General Discussion 

The module STREM2D determines a ‘pseudo’ streamline at the station of interest. Figure 
12 shows a schematic of job stream for STREM2D. The ‘pseudo’ streamline is determined as the 
intersection between the surface geometry and a plane input by the user (PLNST(I), 1=1,4; 
NAMELIST TRAJ). This essentially generates a 2D cut along the surface. This 2D streamline can 
be used for generating data (e.g. pressure coefficient, collection efficiency, heat transfer coefficient) 
for swept and unswept comparisons or for evaluating the traditional quasi-3D icing calculation (i.e. 
calculating swept 3D cases by using 2D calculations along planes normal to the leading edge). This 
module is optional (IRUN=1,5,7,9: and ISTRF=1: NAMELIST TRAJ). If module BSTREM or 
LEWICE2D are to be executed then streamline information is required and either STREM3D or 
STREM2D must be executed. 



Figure 12. - STREM2D segmentation tree structure. 


The first step in determining the 2D streamline is to find the lifting strip associated with the 
section of interest As for the 3D streamline case, this strip is the one associated with the closest 
panel to the section of interest. 

The next step is to determine where the specified plane intersects the local lifting strip. The 
points making up the 2D streamline are essentially the points where the panel edge lines (m-lines 
in the strip) intersect the plane. There will be one point on the 2D streamline for every m-line in 
the strip. This number of points in the 2D streamline then corresponds to the number of points in 
the N-lines for the strip or to the number of panels plus one. 

Each point on the 2D streamline is constructed from an m-line and the specified plane. The 
panel edge lines (A, B) are constructed from the 2 comer points of the panel forming the edge line 
(PI, P2). The intersection of the edge line and the plane is calculated in subroutine PINT. This point 
then is the 2D streamline point. This procedure is repeated for every point in the N-line. 
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2. Printed Output 


STREM2D output consists of coordinates, surface distances (measured from stagnation 
zone where positive surface distance denotes lower surface and negative surface distance denotes 
upper surface), pressure coefficients, and surface normals along the streamline and is written to 
units JOBSUM and OUTPUT. 

H. Subroutine BSTREM 

I. General Discussion 

The module BSTREM interpolates the surface collection efficiencies generated by module 
BETAC onto the streamline generated in STREM2D or STREM3D. Figure 13 provides a sche- 
matic of the BSTREM algorithm. This module is optional (IRUN=5,7,9: NAMELIST TRAJ). If 
LEWICE2D is to be executed, then collection efficiencies along the streamlines are required and 
BSTREM must be executed. 



Figure 13. - BSTREM segmentation tree structure. 


Subroutine BSTREM calculates the collection efficiency along the streamline from the sur- 
face impingement data in two different ways depending on the value of NPSEC. If NPSEC = 4 then 
an interpolation for the streamline points is made from the surface collection efficiency informa- 
tion. If NPSEC= 2, then an extrapolation of the surface collection efficiency is made onto the 
streamline points. 

The first step in making the 3D interpolation (NPSEC = 4) is to determine the location of 
the streamline points relative to the matrix of surface collection efficiency points. For each point in 
the surface streamline, a search is made to determine in which collection efficiency cell the point 
lies (fig. 14). If the point does not lie in any of the cells, then the collection efficiency for that point 
is set to zero. This means care must be taken in setting the width of our section of interest. The 
spanwise width of the section of interest must cover the spanwise range of the streamline in its 
entirety within the impingement local impingement limits. 
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Figure 14. - Collection efficiency interpolation onto surface. 

The interpolation procedure used when NPSEC = 4 is basically that described by Kim (ref. 
8) in his 3D trajectory code paper. Given a point on the surface which lies amid the matrix of sur- 
face collection efficiencies, we have the following interpolation scheme for the collection effi- 
ciency at that poin.: 

P(* 5 ) = p(x c (i+ l,j+ 1)) • u v + $(x c (i,j+ 1)) u • (1-v) + ^ 

P (x(i,j)) • (1-k) • (l-v)+pU(i+lJ)) • ( 1 u) • v 


where: 


u=X s »Xj 

v=X s *Xj 


Eq.5 
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The method employed when NPSEC = 2 is an extrapolation technique based on the 
assumption that there is no spanwise variation in collection efficiency. Alternatively, the method 
assumes that surface lines running parallel to the leading edge are lines of constant collection effi- 
ciency. For each point on the streamline, a test is made to determine which two iso-lines of collec- 
tion efficiency the point lies between. If the point is outside of the impingement limits, then the two 
iso-lines will be zero and the collection efficiency will be set to zero. To determine which iso- lines 
the streamline point lies between, we first form a line parallel to the surface which goes through 
the streamline point. This line (A, B) is formed from the slope of the local panel edge (ATRAN) 
and from the point on the streamline (PI). We then loop through the collection efficiency surface 
points searching for the two closest points to this line (P3, P4), and form a line between these two 
points (A2,B2). The minimum distance between the lines and the points where this minimum 
occurs for each of the lines is then calculated (PIN1, PIN2). If the point (PIN1) lies outside of the 
endpoints of the line segment (A2, B2), then the point lies outside of the impingement limits and 
its impingement efficiency is zero. If the point lies within the line segment, then a linear interpola- 
tion of collection efficiency from the collection efficiencies at the segment endpoints is made. This 
collection efficiency is then the collection efficiency for the streamline point. 

2. Printed Output 

BSTREM output consists of streamline coordinates, surface distances (measured from 
stagnation zone where positive surface distance denotes lower surface and negative surface dis- 
tance denotes upper surface), pressure coefficients, surface normals, trajectory tangents and collec- 
tion efficiency information along the streamline and is written to units JOBSUM and OUTPUT. 

I. Subroutine LEWICE2D 

1. General Discussion 

The module LEWICE2D calculates the ice shape at the section of interest using the collec- 
tion efficiency and pressure coefficient information generated above. The methodology of the ice 
accretion is basically that of the 2D LEWICE calculation expanded to three dimensions.Figure 15 
shows a breakdown of the job stream for LEWICE2D.The ice may be grown either in the surface 
normal direction or in the trajectory tangent direction (see subroutine NEWPTS). If LEWICE2D 
is to be executed, then collection efficiency and pressure coefficient information is needed along 
the streamline and hence BETAC, STREM2D or STREM3D, and BSTREM must be executed. 
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Figure 15. - LEWICE2D segmentation tree structure. 
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The ice accretion process consists of the deposition of super-cooled water droplets on an 
aerodynamic surface and the associated mass and energy exchange processes which result in the 
growth of ice on that surface. That process was first modeled by Tribus (ref. 9) and later devel- 
oped into the model currently used in LEWICE, by Messinger (ref. 10). The Messinger model is 
also used in this code and is applied along streamlines, as calculated by the potential flow portion 
of the code. This chapter describes the Messinger model, the application of that model along 
streamlines, and the input and output files used by these subroutines. 


a. Modeling The Ice Growth Process 


The Messinger model of ice growth is based on the premise that all water impinging on 
the surface of interest exchanges energy with the surface and surrounding environment, resulting 
in freezing of some fraction of the incoming water and the remaining water running back along 
the surface. This model of ice growth on a surface exposed to an icing cloud is depicted in figure 
16, showing the relevant mass and energy exchange processes. 


Impinging water 
droplets 



Convection; 

Evaporation 



Unfrozen 

water 


Figure 16. - Ice growth on a surface 


The ice growth is modeled by dividing the surface into control volumes, along streamlines 
determined from the potential flow analysis, and performing a mass and energy balance on each 
control volume. The control volume extends from the ice-water interface out to beyond the 
boundary layer. Evaluation of the control volumes is started at the stagnation point and marches 
along the upper and lower surfaces to the trailing edge. The mass and energy balances at each sta- 
tion are considered quasi- steady processes. The ice growth is thus modeled as a series of steady- 
state processes with the duration of each step and the number of steps determined by the user. 
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The mass balance, depicted in figure 17, is described by the following equation 

m n + m r - m - = m ; Eq. 6 

c 'in e 'out 1 

where m c is mass flow rate of incoming water, m r is the mass flow rate of runback water from 
the previous control volume, m e is the mass flow rate of evaporated water, m r is the mass flow 
rate of water running back to the next control volume, and m i is the mass flow rate of water leav- 
ing the control volume due to freeze-out. 



Figure 17. - Mass balance control volume 


In this mass balance, the incoming water, incoming runback water, and evaporated water flow rate 
are previously calculated quantities. The amount of water changing phase to ice is determined 
from the eneigy balance and any remaining water is considered to move into the next control vol- 
ume. The value of m r is set to zero at the stagnation point, as this is the start of the marching pro- 
cess for both the upper and lower surfaces. 

The energy balance is handled in a similar manner and is depicted in figure 18. 



Figure 18. - Energy balance control volume 
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The equation describing the control volume energy balance is 


”Vw, T + ^rjw, sur (i - 1) 

= ('Vv, + A rjw, sur + "Vi, jur + *V Al + <7* Aj > 


Eq. 7 


where t w T is the stagnation enthalpy of the incoming water droplets, i w sur j) is the enthalpy 
of the water flowing into the control volume from upstream, i v sur is the enthalpy of the vapor 
leaving the control volume due to evaporation, i Wj sur is the enthalpy of the water running back to 
the next control volume, sur is the enthalpy of the ice leaving the control volume, q c is the heat 
transfer due to convection, and q^ is the heat transfer due to conduction. 

The incoming energy due to water droplet impingement and runback are calculated from 
known information. The energy leaving the control volume due to evaporation and convection 
can be calculated independently. The heat transfer due to conduction is not considered in this 
analysis, as the ice layer is considered to act as an insulating surface. This leaves the energy loss 
due to freeze-out and the energy leaving the control volume due to runback as unknowns. In par- 
ticular, the mass flow rates for these two terms are unknown, as was the case for the mass balance. 
This leaves two equations and two unknowns and the system can be solved. The details of the 
evaluation of each of the terms in the energy equation can be found in appendix A of the LEWICE 
Users Manual (ref. 7). A useful concept for evaluation of the nature of the ice accretion being cal- 
culated is the freezing fraction. This is the fraction of the total water coming into the control vol- 
ume that changes phase to ice. The equation defining freezing fraction is 


m r + m r 

c ' in 

This term can also be used to simplify the evaluation of the energy balance. 

The convection heat transfer term plays an important role in the LEWICE3D energy bal- 
ance. It is through this term that the aerodynamics and the roughness levels can influence the de- 
velopment of the ice accretion. Currently, the convection heat transfer is determined from an 
evaluation of the boundary layer growth on the surface, using an integral boundary layer method. 
The pressure distribution determined by the potential flow code is used as input to the boundary 
layer calculation. The boundary layer calculation determines the displacement thickness and the 
momentum thickness. The Reynold’s analogy is used to determine the heat transfer coefficient. 
Roughness is accounted for by a correlation developed by Ruff. (ref. 7) The complete description 
of the integral boundary layer calculation is found in appendix B of the LEWICE Users Manual 
(ref. 7). 

Expanding the terms in the energy equation as described in the LEWICE manual and com- 
bining Eqs.7 and 8 yields the following form of the energy equation 
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m. 


c (T -273.15) + 


~Pw,s X ‘S ' ' 2 

+ m ^«[ C /W<i-i) * T sur(i - 1) ~ 273 - 15 ) ] +^k As 

~ m e ^ c p W ' Sur ^sur ~ 27 3 . 15 ) +^ v ] 

+ [ ( 1 -f) (m c + w r J - m e ] ^ ( T sur - 273.15) 
+/(m c -m riB ) [c p ' s JT sur - 273.15) - Zy] 


+ /t. 


r.V 7 

fp /y» C ^ 

1 sur~ 1 e~ 2c~ 


Pa J 


As 


Eq. 9 


b. Ice Growth Along Streamlines 


The application of this ice growth model in two dimensions is straightforward in the sense 
that the flow direction of runback water is unambiguous. In three dimensions, the flow direction 
of water out of a control volume is not so easily determined. The most rigorous approach would 
be to solve the air- water flow interaction problem, including the possibility of flow over laige 
roughness elements. In this case, considerable computational effort would be required beyond the 
already significant effort of calculating the flow field, droplet trajectories, and ice growth. As 
such, some degree of approximation is appropriate in order to develop a tool which can be useful 
and not require more computational effort than necessary. 


One alternative is to calculate the boundary layer development over the entire three di- 
mensional surface. This approach requires significant computational expenditure. The approach 
taken for this code is to apply a two dimensional strip analysis along streamlines calculated by the 
three-dimensional panel code. There are differences in streamline direction between those deter- 
mined from a boundary layer analysis and those from the panel code. The differences in ice shape 
development caused by use of the panel code results are most likely smaller than the accuracy lev- 
el of the ice growth prediction method, given the geometric resolution limits established by the 
surface grid. 
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Thus, the integral boundary layer calculation is started at the stagnation line, as deter- 
mined from the panel code results. Then the streamline is divided up into control volumes by us- 
ing the intersection of the streamline with the fore and aft panel edges as the boundaries of the 
control volume and a unit length in the spanwise direction as the other dimension of the control 
volume. Then the fi value at the midpoint of the streamline segment is used as the p for the control 
volume. This use of the streamline (3 value brings in the spanwise influence on the particle trajec- 
tory, whereas a simple cut perpendicular to the leading edge would result in a somewhat different 
P distribution for the ice growth calculation. A representative streamline used for an ice growth 
calculation on a swept wing model is shown in figure 19. 



Grid 

Lines 


Streamline 


Figure 19. - Streamline on swept wing surface. 

The control volumes thus correspond to a series of trapezoidal elements stacked side-by- 
side from the leading edge to the trailing edge on both the upper and lower surfaces. Figure 20 il- 
lustrates a series of control volumes for an arbitrary surface, including the streamline path through 
the center of each volume. 

The P value for the element is taken as the P value at the midpoint of the streamline seg- 
ment. The surface area of the bottom face of the control volume is that of a trapezoid (i.e. equal to 
the panel length times the panel width) and thus is equivalent to the corresponding panel area. 
This value is then used to determine the height of the ice accretion, using m ( and the density of ice 
to determine the ice volume, resulting in the following equation for the height of ice deposited. 
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m t -Ar 


d. = — 1 Eq. 10 

ice o A 

”isur 

The thickness is then considered to be uniform over the entire panel for determination of the new 
geometry. The new coordinates for the panel are obtained from the relation, 


= x i + d ice*i 


Eq. 11 


where Xj is the coordinate of the center of the panel in the i-direction and% is the i-component of 
the unit normal vector for the panel. 



Figure 20. - Control volumes for mass and energy balance. 


As the ice thickness increases, there is the possibility that the ice segments will intersect and thus 
this must be accounted for in the determination of the new geometry. Since this is a strip analysis, 
the ice thickness does not vary along the span at a given chordwise location. Therefore, the possi- 
bility of ice growth intersection is limited to the normal and chordwise directions. In that case, the 
line segments corresponding to the top of every other panel are examined for intersection. If the 
intersection is determined to occur, then a new panel is formed with its center halfway between 
the two old panels. This requires determination of the coordinates of the new panel and renumber- 
ing of the panels. This information is then used in subsequent potential flow calculations. 
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2. Printed Output 


Output from LEWICE2D consists of iced streamline coordinates, surface distances, surface 
normals, trajectory tangents, edge velocities, collection efficiencies ice thickness, heat transfer 
coefficient, and static pressure and is written to units OUTPUT, JOBSUM. 

J. Subroutine BODMOD 

1. General Discussion 

The module BODMOD generates the new geometry file (NGEOM) for the iced geometry 
using the ice accretion at each station of interest (generated in LEWICE2D) and the old geometry 
file (OGEOM). Figure 21 shows a schematic of the job stream for BODMOD.The module reads 
the old geometry file N-line by N-line and if no ice accretion is requested for the section (IMO- 
D(ISEC)=0) it is written to the new geometry file (NGEOM). If an ice accretion has been calcu- 
lated for the section (IMOD(ISEC)=ICEC), then the entire section is modified to reflect the ice 
shape calculated for the two closest streamlines and is written into the new geometry file. In 
essence, the old geometry file (OGEOM) is transferred to the new geometry file (NGEOM) and it 
is overlaid only at sections where ice shapes have been calculated. 
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Figure 21. - BODMOD segmentation tree structure. 


Subroutine BODMOD controls the modification of each N-line and is comprised of two 
steps. First it takes the ice thickness distribution from each of the iced streamlines and extrapolates 
or interpolates that distribution upon each N-line and second, it calculates the new N-line using the 
ice thickness distribution and either the surface normals or trajectory tangents. The algorithm uses 
a cubic interpolation of ice thickness as a function of surface distance in the flow direction and a 
linear interpolation of ice thickness as a function of spanwise position in the spanwise direction. 
Once the ice thickness distribution has been determined for each point on the N-line, the new N- 
line is generated by adding the interpolated ice thickness at each point to the old N-line in either 
the surface normal or the trajectory tangent direction as described in the LEWICE2D section. 

The first step in generating the new N-line is to determine the closest iced streamlines to 
the N-line. If only one iced streamline was calculated then the N-line receives the ice thickness dis- 
tribution calculated for that streamline (DICES(I,1) vs. SST(I,1), 1=1 to NPTS). If multiple iced 
streamlines were calculated, then the two closest streamlines to the N-line will be found. This is 
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done by finding the two closest sections of interest to the N-line and using the iced streamlines 
associated with these sections. The two closest streamlines are denoted by IST1, IST2 and the ice 
thickness distribution and surface distance associated with each of the streamlines are stored in the 
arrays (DICE(I,IST1), SST(I,IST1), I = 1 to NPTS1, DICES(I,IST2), SST(I,IST2), 1=1, NPTS2), 
where NPTS1, NPTS2 arc the number of points in each of the iced streamlines, respectively. In 
addition to the above parameters the parametric distance of the N-line along the line connecting the 
two closest sections of interest is calculated (TL). This distance is used in the linear spanwise inter- 
polation or extrapolation done later. 

The second step is to normalize all of the streamline surface distance arrays to the surface 
distance arrays for the current N-line. This is done to avoid problems in the cubic spline interpola- 
tion which follows. 

The third step involves interpolating the ice thickness distribution from each of the iced 
streamlines onto the N-line surface distance distribution. The interpolation is done using a cubic 
spline fit of the ice thickness vs. surface distance array. 

The fourth step is to determine the ice distribution at the current N-line from the associated 
iced streamlines. If only one streamline is calculated then the N-Line will receive the ice thickness 
distribution from that streamline. If more than one streamline is calculated then a spanwise linear 
interpolation will be made from the two closest ice streamlines (ISTR1, ISTR2). This linear inter- 
polation can be written 


DICE (/) = (DICES (IST2,I) -DICES {IST \, /)) TL + DICES (ISTl, I) Eq. 12 


where TL is the parameteric distance of the N-line along the line connecting the two closest sec- 
tions of interest (ISTR1, ISTR2), and DICE is the ice thickness distribution for the current N-line. 

The fifth and final step is to generate the new N-line. This is done by adding the ice to the 
old geometry in either the surface normal or trajectory tangent direction. This is the same procedure 
used to generate the iced streamlines, and is discussed in the LEWICE2D section. 

2. Printed Output 

Output for BODMOD is written to unit NGEOM and contains the new iced geometry in 
DUGLIFT format. For an explanation of the information in the NGEOM file see table I. 

K. Summary Of Subroutines 

The following tables describe the subroutines used to do the flow, trajectory, streamline, 
impingement efficiency, impingement efficiency interpolation, ice accretion, and geometry modi- 
fication. 
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TABLE L - FLOW FIELD CALCULATION SUBROUTINES 





AIJMX 

SIGMAL 

Computes the matrix, Aij, of dot products of source induced 
velocities with normal vectors on the on-body elements. 

BVORTX 

CONTRL 

Calls PKUTTA or FKUTTA to calculate vortex strength per unit 
path length around the k* lifting strip, B®, for all lifting strips. 

CKARRY 

CONTRL 

Cross checks storage array capacities. 

COLSOL 

SIGMAL 

Solves linear equation matrices for element source densities. 

CONTRL 

DUGLFT 

Controls flow of the modified Hess code which processes body 
surface data for use by the flow and trajectory codes. 

DATPRS 

INPUTL 

Translates, scales, and rotates about the y-axis surface descrip- 
tion data immediately after input 

DKEKFK 

PISWIS 

Calculates D*, E^ and Fufor use in calculating piecewise linear 
spanwise variation of B^*\ 

FAR3 

FLOVEL 

Vectorized subroutine which calculates the induced velocity 
from a lifting panel in the far field. 

FARNL 

FLOVEL 

Vectorized subroutine which calculates the induced velocity 
from a non-lifting panel in the far field.. 

FKUTTA 

BVORTX 

Computes vortex strength, B^, for each lifting strip by the flow 
tangency method. 

FLOVE2 

TRAJEC 

CONFAC 

ARYTRJ 

FLOPNT 

Returns flow velocity for a given point in space. 

FLOVEL 

TRAJEC 

CONFAC 

ARYTRJ 

FLOPNT 

Returns flow velocity for a given point in space (vectorized ver- 
sion). 

HEADER 

CONTRL 

INPUTL 

Writes a printout header. 
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TABLE I. - Continued. 


Subroutine Called bv Description 



NOLIFT 

LIFT 

VIJMX 

PNTVIJ 

DKEKFK 

UNIFLO 

SIGMAL 

AIJMX 

NIKMX 

VELOCY 

PRINTL 


INPUTL 

CONTRL 

Inputs surface quadrilateral comer coordinates and controls 
computation of the geometric properties of lifting quadrilateral 
elements. Also prints the first major DUGLFT output. 

LIFT 

INPUTL 

Computes geometric properties of lifting quadrilateral elements. 

MIS1 

PKUTTA 

FKUTTA 

Linear equation solver. Used in calculation of vortex strengths, 
B®, of lifting strips. 

NEAR 

VFMLFT 

Computes source and vortex induced velocities from an element 
in a lifting strip using the near-field equations. 

NEAR5 

FLOVEL 

Vectorized subroutine which calculates the induced velocity 
from a lifting panel in the near field. 

NEARF 

VFLCFT 

Computes source and vortex induced velocities from an element 
in a lifting strip using the near-field equations. 

NEARNL 

FLOVEL 

Vectorized subroutine which calculates the induced velocity 
from a non-lifting panel in the near field. 

NIKMX 

SIGMAL 

Computes the right hand sides of the Ay matrix, N^ and N 
which are the dot products of the onset flows with the unit nor- 
mal vectors to the on-body elements. 

NOLIFT 

INPUTL 

Calculates geometric quantities for quadrilateral elements in a 
nonlifting section. 
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TABLE L - Continued. 


■JfTTiTWTTmW 

■tnirsrr 


PATPRS 

PINPUT 

Translates, scales, and rotates about the y-axis surface descrip- 
tion data immediately after input. 

PEADER 

PINPUT 

Writes a printout header. 

PISWIS 

VMATRX 

Calculates Vj^ andT^ 1 ^ according to equation (11.2) of ref- 
erence2, and calls DKEDKFK and PSONST to calculate vortex 
induced onset flows. Used only for piecewise linear spanwise 
variation of vortex strength. 

PKUTTA 

BVORTX 

Computes vortex strength, B^, for each lifting strip by the pres- 
sure equality method. 

PNTVIJ 

VMATRX 

If so requested (for debugging purposes only), prints all source 
induced velocities,^;;, and all vortex induced velocities, Vj:^ 
andVjj®. J 

PRINTL 

VELOCY 

Prints the final output of the DUGLFT computations. 

PSONST 

PISWIS 

Computes vortex induced onset flows when the piecewise-linear 
method of spanwise variation of vortex strength is used. 

READ1 

SIGMAL 

COLSOL 

PKUTTA 

FKUTTA 

SUMSIG 

VELOCY 

Reads in one singly subscripted array from a peripheral storage 
unit. 

READ 3 

PNTVIJ 

STEPFN 

PISWIS 

PSONST 

UNIFLO 

AIJMX 

NIKMX 

Reads three subscripted arrays from a peripheral storage unit. 

SETFLO 

FLOPNT 

ARYTRJ 

CONFAC 

TANTRA 

Reads DUGLFT output data stored on unit 18 that is required by 
SR FLOVEL for velocity calculations. If flow velocities are cal- 
culated by other than die Hess method, this code must be 
replaced with a dummy call. 
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TABLE I. - Continued. 


Subroutine Called-bx 


Description 


SIGMAL 

CONTRL 

Controls calculation of the element source densities, cP^ and 
a/~>- 

Computes vortex induced onset flows when the step function 
method of spanwise variation of vortex strength is used. 

STEPFN 

VMATRX 

STOR18 

CKARRY 

Store control and quadrilateral geometrical property data on 
storage unit 18 for use by the flow and trajectory codes. 

SUMSIG 

CONTRL 

Computes the combined element source strengths, Gq-). 

UNIFIX) 

VMATRX 

Stores uniform onset flow velocities for use in calculating ele- 
ment source densities. 

VELOCY 

CONTRL 

Computes the final velocity at the centroid of each element, and 
controls the final printout of the DUGLFT calculation. 

VFLIF3 

FLOVEL 

Vectorized subroutine which calculates the induced velocity 
from a wake panel. 

VFLIFT 

FLOVE2 

Controls computation of velocities induced at a point in space by 
elements of unit source density and unit vortex strength in a lift- 
ing section. 

VFMLFT 

VUMX 

Controls computation of velocities induced at all control points 
by elements of unit source density and unit vortex strength in a 
lifting section. 

VFMNLF 

VUMX 

Computes velocities induced at all control points by elements of 
unit source density in a nonlifting section. 

VFNLFT 

FLOVE2 

Computes velocities induced at a point in space by elements of 
unit source density in a nonlifting section. 

VIJMX 

VMATRX 

Controls computation of source induced velocities,^, and vor- 
tex induced velocities, and at the centroids of all ele- 

ments. 

VMATRX 

CONTRL 

Subexecutive code for computation of induced and onset flow. 


53 



TABLE IL - STREAMLINE CALCULATION SUBROUTINES 



SKtnTrfTfil 


CLINE 

STAGF 

Produces a set of equi-spaced points between the two endpoints 
of a line segment. 

CROSP 

STAGF 

Determines the cross product of two vectors. 

DISLN 

INSTRM 

Computes the minimum distance between two lines and the 
points on each of the lines where this minimum occurs. 

DSTPLN 

INSTRM 

AREAP 

Calculates the minimum distance between a point and a line and 
the point along the line where this minimum occurs. 

FLOVEL 

INSTRM 

STAGF 

STINT 

STREM2D 

STREM3D 

STRMLN 

Vectorized velocity subroutine. Returns a velocity at a given 
point. Description of dependent subroutines given in Tablel. 

FLOVE2 

INSTRM 

STAGF 

STINT 

STREM2D 

STREM3D 

STRMLN 

Returns a velocity at a given point. Description of dependent 
subroutines given in Table I. 

INSTRM 

STREM3D 

Projects off-body integrated streamline points to on-body 
streamline points. 

PANMIN 

STREM2D 

STREM3D 

Determines the closest panel to a given point 

PLIN 

INSTRM 

STREM2D 

AREAP 

Given three points, determines the parametric equation of a line 
parallel to two of the points and which passes through a third 
point. 

STAGF 

STREM3D 

Determines the off body stagnation point at a section of interest. 

STINT 

STRMLN 

4th order Runge-Kutta integration scheme for calculating 
streamlines. 
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TABLE II. - Concluded. 


Subroutine Called by Description 


STREM2D 

MAIN 

Calculates a 2D streamline based on the intersection of a given 
plane and the geometry. 

STREM3D 

MAIN 

Controls the calculation of the 3D integrated on-body stream- 
lines. 

STRMLN 

STREM3D 

Integrates a 3d streamline from a given point 

TRANSF 

INSTRM 

STREM2D 

Projects a set of points onto a given plane along a specified direc 
tion. 
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TABLE ffl. - COLLECTION EFFICIENCY CALCULATION SUBROUTINES 




Description 

ACORR 

BETAC 

Determines surface area correction due to sweep of surface for 
the “quasi” 2D impingement efficiency calculation. 

ARYTRJ 

BETAC 

Calculates trajectories for a given set of release points. 

AREAP 

BETAC 

Calculates the area of an N-sided polygon. 

BETAC 

LEWICE3D 

Controls the calculation of surface impingement efficiency at a 
section of interest 

CDRR 

PRFUN 

PARTCL 

Given Reynolds number, returns Davies number for a sphere. 
Used for water drops for which Reynolds number is less than 
equal to 81.23. 

CLINE 

BETAC 

Produces a set of equi-spaced points along a line between two 
given points. 

CONFAC 

BETAC 

Determines the release point for a trajectory which passes within 
a specified tolerance of a target point at the section of interest. 

DVDQ 

TRAJEC 

Integrates particle equations of motion for each time step. 

FALWAT 

PARTCL 

Returns still-air, terminal settling speed for a water drop. Uses 
equation of Beard. 

IMPACT 

TRAJEC 

Used in runs under control of CONFAC to adjust trajectory ini- 
tial y, z coordinates to avoid impact on the body on the next tra- 
jectory after impaction has occurred. This is a problem-specific 
subroutine that must be programmed by the user. 

IMPLIM 

BETAC 

Determines tangent trajectories at the section of interest. 

MAP 

CONFAC 

Controls the iterative calculation of trajectories to a specified tar- 
get point. 

MATINV 

MAP 

Linear equation solver. 

PARTCL 

BETAC 

Reads particle specification data and returns still-air, terminal 
particle settling speed and other particle data as required for the 
particular type of particle. This is a particle type-specific code. 
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TABLE in. - Continued. 







The version provided here is for water drops. 

POLYGO 

CONFAC 

Calculates area of plane polygon of N vertices. Provides cross- 
sectional areas of particle flux tubes which are used to compute 
concentration factors, concentration ratios and collection effi- 
ciencies. 

PRFUN 

TRAJEC 

Given the particle Reynolds number, returns the factor which 
when multiplied by^Vp - V a yields the first term on the right side 
of eq. (1). This is a particle type-specific function. The version 
provided here is for water drops. 

STRPNT 

TANTRA 

Specifies a curve in three-dimensional space on which lie the ini- 
tial points of all trajectories used in computing a tangent trajec- 
tory to the body. Also specifies coarse and fine step sizes to be 
used in traversing the curve in search of the tangent trajectory, 
and it steps along the curve to define new initial trajectory points 
under control of TANTRA. The version supplied here uses 
straight line curves. 

TRAJEC 

ARYTRJ 

IMPLIM 

MAP 

Computes particle trajectories. 

TRANSF 

PRFUN 

MAP 

Transforms coordinate system for the “flow system” to the “flux 
tube system”, or reverse. 

WCDRR 

PRFUN 

PARTCL 

Given Reynolds number, returns Davies number for a water 
drop. Used for cases where the Reynolds number is greater than 
81.32. 

BETINT 

BSTREM 

Determines which surface collection efficiency cell the stream- 
line point lies in if any and the weighting factors for the interpo- 
lation. 

BSTREM 

LEWICE3D 

Controls the interpolation of surface collection efficiency onto 
the streamline. 

DISLN 

BSTREM 

Determines the minimum distance between a point and a line 
and the point on the line where this minimum occurs. 
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TABLE m. - Concluded. 



ImiiiimiliR 


.umiraiir 




PLIN BSTREM Given three points determines the parametric equation of a line 

parallel to two of the points and which passes through the third 
point. 
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TABLE IV. - ICE ACCRETION CALCULATION SUBROUTINES 


Subroutine Called bv Description 


BDYLR 

ICECAL 

Determines heat transfer coefficients and transition points for 
streamline. 

COMPF 

EBAL 

Solves the energy equation for the freezing fraction. 

COMPS 

LEWICE2D 

SEGSEC 

Calculates surface distance for streamline. 

COMPT 

EBAL 

Solves the energy equation for the ice surface temperature given 
a value of freezing fraction. 

CONST 

LEWICE2D 

Sets constants in /ICECOM/ for ice accretion calculation. 

CPW 

COMPF 

COMPT 

Calculates specific heat of water for a given temperature. 

DSTPLN 

BDYLR 

Determines the minimum distance between a point and a line 


COMPS 

ICECAL 

NWFOIL 

NWPTS 

and where this point on the line occurs. 

EBAL 

ICECAL 

Controls the mass and eneigy balance for each of the segments 
on the streamline. 

ICECAL 

LEWICE2D 

Controlling routine for ice distribution thickness and new airfoil 
point calculations at each step. 

INTRST 

SEGSEC 

Determines if two line segments in a line intersect and if so, at 
which point this intersection occurs. 

LEWICE2D 

LEWICE3D 

Controls the ice accretion calculation for a streamline. 

NWFOIL 

ICECAL 

Computes the new x, y, z coordinates for the iced airfoil in the 
surface normal or trajectory tangent direction. 

NWPTS 

LEWICE2D 

Tests the iced streamline point distribution for refinement. If seg- 
ments have become too large dining a time step, they are subdi- 
vided into two segments of equal size. 
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TABLE IV - Concluded. 





PLIN 

INTRST 

Given three points determines the parametric equation for a line 
parallel to two of the points and which passes through a third 
point. 

PLNFRM 

SEGSEC 

Produces a plane given three points. 

PLNNRM 

NWFOBL 

Calculates cross product of two vectors. 

PVI 

COMPT 

Calculates the vapor pressure over ice for a given temperature. 

PVW 

COMPF 

COMPT 

Calculates vapor pressure over water for a given temperature. 

RHOICE 

ICECAL 

Calculates local ice density using Macklin correlation (ref. 7). 

SEGSEC 

ICECAL 

Removes segments that intersect due to ice growth. 

TRANSF 

NWFOIL 

SEGSEC 

Projects a set of points onto a given plane along a specified direc- 
tion. 

VEDGE 

LEWICE2D 

Determines edge velocity, static temperature, and pressure along 
a streamline. 
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TABLE V. - GEOMETRY MODIFICATION CALCULATION SUBROUTINES 




Descriotion 

BODMOD 

LEWICE3D 

Controls the geometry modification at each ice step. 

CSPLINE 

NLNMOD 

Cubic spline interpolation routine. 

DSTPLN 

ISCFND 

Determines the minimum distance between a point and a line 
and where this point occurs on the line. 

GEOMOD 

BODMOD 

Controls the modification of each N-line on a lifting section and 
its subsequent loading into the new geometry file. 

ISCFND 

NLNMOD 

Finds the two closest iced streamlines to a given streamline and 
its relative position between them. 

NLNDAT 

NLNMOD 

Determines the surface normals for each point on a given N-line. 

NLNMOD 

GEOMOD 

Controls the calculation of the ice thickness distribution and 
geometry modification for each N-line. 

NORM 

NLNMOD 

Normalizes an array of surface distances to a given surface dis- 
tance. 

NWFOI2 

NLNMOD 

Calculates the new points on the N-line given the old N-line 
points, ice thickness distribution and either the surface normals 
or the trajectory tangents. 

PANMIN 

BODMOD 

Determines the number of the panel closest to a given point. 

PLIN 

ISCFND 

NLNDAT 

Given three points, determines the parametric equation of a line 
which is parallel to two of the points and which passes through 
the third. 

PLNFRM 

GEOMOD 

NLNMOD 

Forms a plane given three points on the plane. 

SURFD 

NLNMOD 

Determines the surface distance distribution from a set of points. 

SWITC1 

NLNMOD 

Transfers points from a two-dimensional array into a one-dimen- 
sional array. 

SWITC2 

NLNMOD 

Transfers points from a three-dimensional array into a one- 
dimensional array. 
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TABLE V. - Concluded. 


Subroutine Called bv 


. P esc r ip-tipn 


TRANSF GEOMOD 

NLNMOD 

TRIB CSPLINE 

WEIGT NLNMOD 


Projects a set of points onto a given plane along a specified direc- 
tion. 

Solves for coefficients in a cubic spline curve fit. 

Determines the ice distribution at an N-line from weighted val- 
ues of the two closest iced streamlines. 
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III. INPUT FILES 


Two basic input files are required to run the code and a third is optional if the restart capa- 
bility of the code is used (IRES=1). A geometry file (unit OGEOM) is required for the flow field 
generation. The job control file (unit INPUT) is required and contains flags and inputs for the tra- 
jectory, and ice accretion calculations. The third file, which is optional, is a restart file (unit 
RESTRT) which allows the user to continue from the point where the last run was terminated. This 
file is useful for long runs where it might be more advantageous to split the job into smaller runs. 
A brief description of the flow field input file is contained in Hillyer Norment’s trajectory code 
manual (ref. 5) while a more detailed description is available in the Duglift users manual (ref. 4). 

The run parameter file (unit INPUT) contains basically three namelists which control the 
trajectory and ice accretion calculations. A description of each of the variables and namelists is 
given in table VI. 

The DUGLIFT flow field input file (unit NGEOM) contains geometry information in DUG- 
LIFT format. Table VH, which was taken from reference 5, gives a brief description of the input 
format of the variables. 

In addition to the flow field and job control input files there is an optional restart file. This 
file allows a job to be restarted from its previous termination point. To restart a job the restart flag 
must be set (IRES=1) and the previous restart file must be provided on unit REST. This file contains 
collection efficiency and ice accretion information at each time step and section of interest. The 
restart file is read in subroutine REST. For information about the type and format of the data see 
subroutine REST. 
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TABLE VI. - LEWICE3D STANDARD INPUT FILE DESCRIPTION. 
NAMELIST Variables an d Format Description 


IMPING IRUN, IFLOW, ICE, IRUN 
ISTRF, ICEC, IRES, 

NPSEC, NBR, NBC, 

XSEC, YSEC, ZSEC, =1 

xsa, Ysa, zsa, 

XTIP, YTIP, ZTTP, 

DSHIFT, SHFTF, 

PLNST (NAMELIST =2 

FORMAT) 

=3 


=4 


=5 

=6 


=7 


=8 


Rag controlling trajectory and streamline 
calculations. 

Only streamline calculation will be carried 
out. Must input XSEC, YSEC,ZSEC. (Sub- 
routines STREM2D or STREM3D will be 
called depending on flag ISTRF). 

CONFAC run will be carried out to deter- 
mine trajectories that pass XSEC, YSEC, - 
ZSEC, (subroutine CONFAC). 

Tangent trajectories will be determined. 
Must input XSCI,YSCI,ZSCI. (subroutine 
TANTRA). 

An array of particles will be released and col- 
lection efficiencies will be calculated. Must 
input XTIP, YTIP, ZTIP. (subroutine 
ARYTRJ). 

All of the above will be calculated. Must 
input XSEC, Y SEC,ZSEC. 

CONFAC run will be carried out, followed 
by a tangent trajectory, and a collection effi- 
ciency run. 

Streamline calculation will be carried out fol- 
lowed by a tangent trajectory, and a collection 
efficiency run. Must input XSEC, YSEC, 
ZSEC. 

Tangent trajectory calculation will be fol- 
lowed by a collection efficiency run. Must 
input XSCI, YSCI, ZSCI. Collection effi- 
ciency run will be carried out and will be fol- 
lowed by tangent trajectory run. Must input 
XSEC, YSEC, ZSEC. 


EFLOW Row field control flag. 

=0 Row solver will not be run. Must provide 
geometry on OGEOM. 

=1 Row solver will be run. Must provide geom 
etry on unit OGEOM. 

ICE Ice accretion calculation control variable. 
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TABLE VI. - Continued. 
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=0 Lewice2D ice accretion calculation will not 
be run. 

=1 Lewice2D ice accretion calculation will be 
run. Must provide accretion calculation vari- 
ables. (NAMELIST ICEIN) 

ISTRF Streamline calculation control variable. 

=0 A 3D streamline will be integrated along the 
surface at the section of interest. Must input 
XSEC, YSEC, ZSEC. 

=1 A 2D cut will be generated along the surface. 
The 2D slice will be the intersection between 
the surface and the plane input by the user 
(PLNST(ICEC,4» where the plane is: 

PLNST(ICEC, 1 )*X + PLNST(ICEC,2)*Y + 
PLNST(ICEC,3)*Z + PLNST(ICEC,4) = 0 

The user must input PLNST, and XSEC, 
YSEC, and ZSEC. 

ICEC The number of sections for which the above 
trajectory or ice accretion calculations will be 
made. 

IRES Restart flag. 

=0 No restart will be made and job will run from 
the beginning. 

=1 Job will continue from last point of execu- 
tion. Must link restart file (unit 26) 

NPSEC Variable controlling the type of region at the 
section of interest 

=2 The region at the section of interest is a line 
and hence only two points are needed to 
describe it (i.e. XSEC(ICEC,1), YSEC(I- 
CEC,1), ZSEC(ICEC,1 ) and XSEC(ICEC,2), 
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TABLE VL - Continued. 


NAMELIST Variables and Format Description : 

YSEC(ICEC,2), ZSEC(ICEC,2) fully 
describe the region at the surface). This type 
of calculation is justified for regions where no 
spanwise variation in the flow field or collec- 
tion efficiency is expected. A single row of 
trajectories will be released along the section 
line and a 2D beta calculation will be used for 
determining collection efficiency. Flow field 
data are linearly extrapolated onto the stream- 
line assuming no spanwise variation. 

=4 The region at the section of interest is a rect- 
angle and hence four points are needed to 
describe the region of interest (i.e. XSEC(I- 
CEC, 1), YSEC(ICEC, 1), ZSEC(ICEC, 1) and 
XSEC(ICEC,2), YSEC(ICEC,2), ZSEC(I- 
CEC,2), XSEC(ICEC,3), YSEC(ICEC,3), 
ZSEC(ICEC,3) and XSEC(ICEC,4), 
YSEC(ICEC,4), ZSEC(ICEC,4) describe the 
four comers of the rectangle at the section of 
interest. This type of calculation is for regions 
where the flow is expected to be fully 3D. A 
matrix of trajectories will be released into the 
rectangle of interest to generate a distribution 
of collection efficiencies on the surface. A 3D 
collection efficiency is made. The collection 
efficiency and flow field data are interpolated 
onto the streamline using linear interpolation. 

NBR The number of rows of trajectories to be 

released at each section of interest 
NBR(ICEC). Typical value is 20. 

NBC The number of columns of trajectories to be 
released at each section of interest NBC(I- 
CEC). For the 2D Approximation (i.e. 
NPSEC(ICEC)= 1 ) NBC will be set to one 
and only a line of NBR(ICEC) trajectories 
will be released at the section of interest. 


XSEC, 
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TABLE VI. - Continued. 


NAMELIST Variables and Format Description — 

YSEC, 

ZSEC Arrays describing the region of interest. 

Depending on NPSEC(I), 1=1, ICEC either a 
line is desired (NPSEC=2) and two points 
along this line must be entered, or a rectangle 
is desired (NPSEC=4) and the four comer 
points of the rectangle of interest must be 
entered. The points must be off-body points. 
These arrays are needed to run the streamline 
and CONFAC calculations. 

XSCI, 

YSCI, 

ZSCI Arrays either generated by subroutine CON- 
FAC or input by the user that define upstream 
release points for trajectories that pass 
through the points XSEC, YSEC, ZSEC at 
the region of interest These arrays are needed 
to run the tangent trajectory routine. 

XTIP, 

YTIP, 

ZTIP Arrays either generated by subroutine 

TANTRA or input by the user that define 
upstream release points for tangent trajecto- 
ries for the upper and lower surface along the 
line defined by XSCI, YSCI, ZSCI at the 
region of interest. These arrays define the 
region to release impacting trajectories. 
These arrays are needed to run the ARYTRJ 
trajectory subroutine which generates collec- 
tion efficiency data. 

DSHIFT Normal distance off-body where the stream- 
line integration is started. Because of velocity 
gradient problems at panel edges, integration 
of a streamline at the surface is difficult. For 
this reason the streamline integration is made 
at a distance off the body equal to DSHIFT. 
Typical values are 0.002. 
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TABLE VL - Continued. 


NAMELIST Variables and Format Pescriptfoa 


TRAJ 


SHFTF Variable which controls the amount the sur- 
face is shifted to overcome difficulties in inte- 
grating the trajectories due to high velocity 
gradients near panel edges. The surface is 
shifted in the flow direction an amount equal 
to SHFTF*DHSEFT. The default value for 
SHFTF is 0.0. Typical values may range from 
0.0 to 1.0. 

PLNST Array defining plane which is to cut surface 
to generate 2D streamline at each section of 
interest. Plane is defined as 


PLNST(ICEC, 1 )*X + PLNST(ICEC,2)*Y + 
PLNST(ICEC,3)*Z - PLNST(ICEC,4) = 0 

Array PLNST must be entered if ISTRF=1. 

IPLOT Logical variable controlling output of trajec- 
tory information. 

IPLOT, VINF, CHORD, 

RHO, VIS, HI, HMINI, =TRUE Trajectories are written to unit TEMP24. 
IDIS,TPRINT, 

XSTART, XFINAL =FALSE No trajectory data is output. 

EPSE, NW, RW, TOL, 

XE, YE, XI, YI,DMDS, VINF Airspeed (M/S). 

PLWC, FNR, DFINE, 

(NAMELIST FOR- CHORD Characteristic dimension of body (M). 

MAT) 

RHO Air density (Kg/M3). 


VIS Air viscosity (Kg/(M-S). 

HI Initial time step for numerical integrator. 

Typical value is 0.1. 

HMINI Minimum time step for numerical integrator. 
Typical value is 0.00001. 
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TABLE VL - Continued. 
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IDIS The number of the particle sizes in the parti- 
cle distribution. If IDIS is greater than one, a 
particle distribution is assumed and 
DMDS(I), PLWC(I), 1 = 1, IDIS must be 
input. If IDIS = 1 than PWLC(1)=1.TPRINT 
Output time interval for trajectory plotting 
arrays. Typical value is 0. 1 . 

XSTART Initial X release plane for trajectory calcula- 
tions (non-dimensional). 

XFINAL Termination X plane for trajectory integra- 
tion (non-dimensional). 

EPSI Array used to control local error in trajectory 
integration in each of the coordinated direc- 
tions. Typical values are 0.000001. 

NW Number of trajectories used to define the flux 

tube periphery. This parameter should only 
be greater than one if off body concentration 
factors are desired. If NW =1, then single tra- 
jectories are computed. 

RW Radius of particle flux tube in target plane. 
Only used if NW is greater than one. 

TOL Tolerance for reaching a point on tangent 

plane. Controls how closely trajectories pass 
through points XSEC, YSEC, ZSEC in the 
CONFAC calculation. 

XE, YE Target point coordinates of the last three 
guesses. Used in subroutine CONFAC in 
search for target point trajectories (see sub- 
routine CONFAC). 

XI, YI Initial point coordinates of the last three 

guesses. Used in subroutine CONFAC in 
search for target point trajectories (see sub- 
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TABLE VL - Concluded. 


NAMELIST Variables and Format Description 


DMDS 


PLWC 


FNR 


DFINE 

LWC 

ICEIN LWC, TAMB, PAMB, TAMB 

RH, XKINIT, SEGTOL, 
QCOND,TSTOP, PAMB 

DTIME, DTFLW 
(NAMELIST FOR- RH 

MAT) 

XKINIT 

SEGTOL 


QCOND 

TSTOP 

DTIME 


DTFLW 


routine CONFAC). 

Distribution of droplet or ice aggregate diam- 
eters to be run: DDS(I), 1=1, IDIS. 

Distribution of percent liquid water content 
for the distribution of particles. If IDIS =1 
then PLWC = 1. 

Froude number. If gravitational forces in the 
z-direction are to be considered then FNR = 
1. The default value for FNR is 0 or no grav- 
itational forces. 

Step size used in search for tangent trajectory. 

Liquid water content of cloud (g/m ). 

Ambient temperature (K). 

Ambient pressure (Pa). 

Relative humidity of cloud (percent). 

Roughness factor (m). 

Maximum growth length of a surface seg- 
ment before it is divided into two surface ele- 
ments. 

Length of icing encounter (sec). 

Time stepping for ice growth (sec). Should be 
fraction of DTFLW. 

Time stepping for flow field calculation (sec). 
Should be fraction of TSTOP. 
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TABLE VIL - DUGLIFT INPUT FILE DESCRIPTION 


Card No. Variables and Format Description 


Card No. Variables and format Description. 


1 

2 


TITLE(I), 1=1, 18), A4 Run identification. 

CASE, LIFSEC, Run control data: 

IATACK, NSYM1, 


NSYM2, MPR, LEAK, CASE 

FRAC, MACH (col. 1-4) 

(A4,6I4,2X,2F 10.0) 

LIFSEC 
(col. 5-8) 

IATACK 
(co. 9-12) 


NSYM1 
(col. 13- 16) 
NSYM2 
(col. 17-20) 


Four character body identification. 


Total number of lifting sections. 

Number of angles of attack (i.e., uniform 
free stream flows) to be specified via cards 
no. 4. Maximum value is 10. If the com- 
pression correction is to be applied (MACH 
> 0.0), it is necessary that IATACK = 1. 

One of the three values 0, 1 or -1 is entered. 
NSYM1 specifies the 1st symmetry plane, 
and NSYM2 specifies the second 
symmetry plane according to 


0 nonexistent. 

+1 a plus (ordinary symmetry plane. 

- 1 a minus (anti) symmetry plane. 

MPR Print flag used for program debugging only, 

(col. 21-24) 

0 No debug print. This is the normal 
value for this parameter. 

1 Print the source induced velocity 
matrix, Vjj, and, if LIFSEC > 0, 
print the dipole induced velocity 
matrices, and Vy^. 

>2 Print the dot product matrices Ay, 
Nj^ and n/°°\ and the element 
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TABLE VII. - Continued. 


Card No. Variables and Format 


3 IPROS, LOFF, 

MOMENT, LIST, IOUT 
(5LI) 


Description 

source densities, of® and of 00 ' 

3 Print the onset flow matrices, 
and Vj^ 


LEAK Number of inlet quadrilateral elements, 
(col. 25-28) These must be the first elements in the dig- 
ital description set (cards no. 12). 

FRAC Fraction of unit free stream flow that passes 
(col.3 1 -40) through each of the LEAK inlet elements. If 
LEAK = 0, leave this field blank. 


MACH Mach number of the free stream flow. (Note 
(col.41-50) that this is a floating point number.) If N^ 
< 0.5, leave this field blank. 


Logical control flags: 

IPROS If true, the card 12, 13, and 15 coordinates 
(col. 1) are to be translated, scaled, and rotated 

about the y axis before processing, and card 
no. 5 is to be input. If false, no translation, 
scaling or rotation is done, and card no. 5 is 
not input. 

LOFF if true, velocities at off-body points are to 
(col. 2) be calculated. The off-body points are spec- 
ified by the user via input of the no. 15 
cards. If false, off-body velocities are not 
calculated and there is no input of the no. 1 5 
cards. 

MOMENT If true, the moment origin is specified by 
(col.3) input of card no. 6. If false, card 6 is not 
input and moments are computed about 
point (0,0,0). 

LIST If false, specifies complete execution. If 

(col. 4) true execution is terminated after the first 
main part of the printed output. 
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TABLE VII. - Continued. 


Card No. Variables and Format Description 


4 


5 


6 

7 


(ALPHAX(I), 
ALPHAY(I), ALP- 
HAZ(I), 1=1, IATACK) 
(3E10.0) 


IOUT If true, the 29 geometric quantities for each 
(col. 5) nonlifting element and 45 geometric quan- 
tities for each lifting element are printed. 
The normal value for this parameter is 
false. 

Direction cosines of uniform onset (i.e., free stream) 
flow vectors. IATACK is the number of uniform onset 
flows specified in card no. 2. One set of direction cosines 
per card. If the compression correction is applied 
(MACH . 0.0), only one uniform onset flow vector can be 
specified. If more than one vector is specified, only the 
first is passed along via unit 18 for use by SR SETFLO 
and FLOVEL. The direction cosines are with reference 
to the airplane coordinate system (after rotation by angle 
ANGLE (card 5)). These vectors are equivalent to unit 
free stream velocities. Ordinarily, free stream unit veloc- 
ity components are (1.0, 0.0, 0.0). 


ANGLE, XSCALE, 

YSCALE, ZSCALE, Input only if IPROS = TRUE on card 3. Same as card no. 

XTRANS, YTRANS, 4 of program PBOXC. 

ZTRANS, (7F10.0) 


ORIGNX, ORIGNY, 

ORIGNZ (3E10.0) Coordinates of the moment origin. This card is input only 

if MOMENT = TRUE on card 3. 

LKUTT, LASWAK, 

PESWIS, IGW (5L1) Input only if LIFSEC > 0 on card 2. Logical control flags 

for lifting section data 


LKUTT If true, the flow tangency method for 
(col. 1) application of the Kutta condition is 

selected. This means that one point in or 
near the wake of each lifting strip (not 
counting extra strips) must be specified via 
input of card no. 9. If false, the pressure 
equality method is selected, and cards no.9, 
13, and 14 are not input 
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TABLE YU. - Continued. 


Card No. Variables and Format 


8 (NSORCE(J), 

NWAKE(J), NSTRIP(J), 
NLINE1(J), NLINEN(J), 
IXFLAG(J), J=l, LIF- 
SEC), (614) 


Description 

LASWAK If true, the trailing edge of the last wake 
(col. 2) element is automatically extended by the 
code to x = °°. This is the semi-infinite last 
wake element option. 

PESWIS If true, the piecewise linear method for 
(col. 3) calculating spanwise variation of lift vorac- 
ity is selected, and lifting strip widths must 
be input via cards no. 11. If false, the step 
function option is selected, and cards no. 1 1 
are not input. 

IGW If true, there are ignored lifting elements 

(col. 4) which must be defined via input of the no. 

10 cards. If false, there are not ignored ele- 
ments, and cards no. 10 are not input. 


This card only input if LIFSEC > 0 on card 2. 


NSORCE Number of on-body elements (including 
(col. 1-4) ignored) in each lifting strip of the Jth lift- 
ing section. 

NWAKE Number of wake elements in each lifting 
(col. 5-8) strip of the Jth lifting section, including a 
semi-infinite final element if this option is 
selected. 


NSTRIP Number of lifting strips in the Jth lifting 
(col. 9-12) section.Include extra strips only if they are 
defined via input of cards no. 12. 

NLINE1 If the piecewise linear option is selected, 
(col. 13-16) (PESWIS = TRUE on card 7), NLINEl(J) 
specifies the edge condition of the first strip 
on the Jth lifting section. If the step function 
option is specified, ignore this field. 


NLINEN Same as NLINEl(J) but for the last strip of 
(col. 17-20) the Jth lifting section. 
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TABLE VII. - Continued. 


Card No. Variables and Format Description 

IXFLAG IXFLAG(J) = 0 means that no extra strips 
(col. 21-24) are defined via input (i.e. via cards no. 12). 

IXFLAG(J) = 1 means the first strip is an 
extra strip. If the piecewise linear option is 
selected, this also requires NLINEN(J) = 4 
or 5. 

IXFLAG (J) =3 means the last strip is an 
extra strip. If the piecewise linear option is 
selected, (PESWIS = TRUE on card 7), this 
also requires NLINEN(J) = 4 or 5. 

DCFLAG(J) = 2 means that both the first 
and last strips are extra strips. If the piece- 
wise linear option is specified, this requires 
that both NLINEl(J) and NLINEN(J) = 4 
or 5. 

LIFSEC is specified on card 2. A separate card is 
required for each lifting section, and the cards are input 
in the same order as input of the quadrilateral data via 
cards no. 12. 

Input only if LIFSEC > 0 on card 2 and if LKUTT = true 
on card 7. Number of points which the flow tangency 
method for application of the KUTTA conditions is to be 
applied. It is required that KUTTA equal the total num- 
ber of lifting strips, not counting extra strips. KUTTA is 
used to read the point coordinates, and the unit vectors 
normal to the wake or airfoil surface at these points, via 
cards no. 13 and 14. 

Input only if LIFSEC > 0 on card 2 and if IGW = TRUE 
on card 7.1 = lifting strip index: J = lifting section index. 
If on the Ith strip of the Jth lifting section there is a sub- 
strip of ignored elements, the substrip is defined by spec- 
ifying its beginning and ending element indices via 

IG 1 (I, J) = index of the first ignored element on the 

lifting strip. 


9 KUTTA(I4) 


10 ((IG1(I,J), IGN(IJ), 1=1, 

NSTRIP(J)), J=l, LIF- 
SEC), (1214) 
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TABLE VII. - Continued. 


Card No. Variables and Format Description 


11 


12 


IGN(I,J) = index of the last ignored element on the 
lifting strip. 

If there are no ignored elements on a strip, leave both 
fields blank: but IG1 and IGN must be specified for every 
lifting strip if IGW = TRUE on card 7. Six strips per 
card. Each lifting section begins a new card. 

LIFSEC is specified on card 2, and NSTRIP(J) on card 8. 


(WIDXTR(U), 
(WIDTH(I,J), 1=2, 
NSTRIP(J)-K), WIDX- 
TR(2,J),J=1, LIFSEC), 
(7E10.0) 

K = 0 if IXFLAG( J) = 0 

K = 1 if IXFLAG(J) = 1 or 
3 

K = 2 IF IXFLAG(J) = 2 


Input only if LIFSEC > 0 and if PESWIS = TRUE on 
card 7. These quantities are the widths of each lifting 
strip for use in calculating the spanwise variation of vor- 
ticity via the piecewise linear method. 

WIDXTR( 1 J) specifies the width of the first extra strip 
of the Jth lifting section. If NLINE1(J) 
NE 4, leave this field blank. 

WIDTH(I,J) specifies the width of the Ith lifting strip 
of the Jth lifting section. 

WIDXTR(2J) specifies the width of the last extra strip 
of the Jth lifting section. If NLINEN(J) 
NE 4, leave this field blank. 


LIFSEC is specified on card 2, and NSTRIP(J), 
NLINE1(J), NLINEN(J) and IXFLAG(J) are specified 
on card 8. 


X, Y, Z, STAT, LAB, 
XX, YY, ZZ, STATT, 
LABL, (3E 10.0,212/ 
3E10.0,2I2) 


On-body wake quadrilateral element comer point coordi- 
nates are specified by these cards for both lifting and 
nonlifting sections, one point per card. The body and 
wake surface panels are constructed from these data. 
(Note: there must be an even number of no. 12 cards. 
Add a blank card to the end of the card 12 deck if neces- 
sary). 


X, Y, Z Quadrilateral comer point coordinates. 
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TABLE VII. - Continued. 




jx^oSSjutu 


in 


XX, YY, ZZ 
(col. 1-30) 

STAT Status parameter: Allowed values are 0, 

STATT 1, 2,3: 

(col. 32) 

0 This point is on the same N-line 
as the last point. 


1 This point starts a new N-line. 


2 This point starts a new section. 

3 This is the last point in the card 12 
input. 

LAB Specifies a lifting or nonlifting section 

LABL 

(col. 32) 0 Nonlifting. 


1 Lifting. 

This field is relevant only when STAT 
or STATT = 2, that is, only on the 
first card of a new section. 

13 (XC(I),YC(I),ZC(I),I=1, 

KUTTA) (3E10.0) Input only if LKUTT = TRUE on card 7. KUTTA is spec- 

ified via the card 9 input. Coordinates of points (one point 
per card) at which the Kutta condition is to be applied via 
use of the flow tangency method. If IPROS = TRUE (card 
3), the coordinates according to the card 5 input data. 

14 (XN(I), YN(I), ZN(I), 

1=1 , KUTTA) (3E10.0) Input only if LKUTT = TRUE on card 7. KUTTA is spec- 
ified via the card 9 input. Components of the unit vectors 
(one vector per card) normal to the wake or airfoil surface 
at the points specified by the no. 13 cards at which the 
Kutta condition is to be applied via use of the flow tan- 
gency method. The order of input must be consistent with 
that of the no. 13 cards. If IPROS = TRUE (card 3), a trans- 
formation is automatically applied by the code to adjust for 
rotation of coordinates by angle ANGLE (card 5). 
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TABLE VII. - Concluded. 


Card No. Variables and Format Description 


15 


XOF, YOF, ZOF, STAT, 

XOFF, YOFF, ZOFF, XOF, YOF, ZOF Coordinates of off-body points 

STATT, (3E10.0, 12/ XOFF, YOFF.ZOFF at which flow velocities are to be 

3E 10.0,12) (col. 1-30) calculated, one point per card. 


STAT Status parameter. A value of 3 

STATT signifies the end of the off-body 

(col. 32) points. Otherwise, leave this field 

blank. 


If IPROS=TRUE (card 3), the code automatically trans- 
lates, scales, and rotates these coordinates according to 
the card 5 input data. 
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V. EXAMPLE CASE 


For the test case, a swept NACA 0012 airfoil with a 30 degree sweep, an aspect ratio of 1.4, 
and a .4399 meter chord was used (fig. 22). Calculations at 0 degrees angle of attack were made 
for this configuration at two different airspeeds, 150 MPH and 165 MPH. This choice allowed for 
direct comparison to experimental impingement data for the 165 MPH case and to ice accretion 
data for the 150 MPH case. Because the cases were very similar, only the ice accretion related run 
files have been included. The experimental impingement calculations are summarized in figure 23. 
The various steps in producing an ice shape for a 3D geometry are illustrated. 

The ice accretion calculations were carried out in a single run on the Cray XMP at LEWIS. 
The calculation conditions were airspeed, 150 MPH; angle-of-attack, 0 degrees; drop size, 20 
microns; LWC, .75 g/m 3 ; icing time, 2 minutes. Three sections of interest were chosen to resolve 
the spanwise ice shape. Because of the relative shortness of the icing encounter a single step ice 
accretion calculation was chosen. These sections of interest were located at the 10, 50 and 90% 
span locations. The panel model contained one lifting section containing 14 lifting strips of con- 
stant width and 91 chordwise segments. This model yielded 1218 lifting elements. Figures 24 and 
25 show the job control input file (unit INPUT) and the flow field input file respectively (unit 
OGEOM). Figure 26 contains the job summary file (unit JOBSUM). The entire calculation 
required approximately 1340 seconds on the Cray XMP. 

The first step is to generate a flow field. If the same flow field will be used for several tra- 
jectory runs then it is suggested that the flow field be generated on the first run (IFLOW=l) and 
saved (unit FLOWF) for any runs thereafter (IFLOW=0). It is also a good idea when calculating 
flow for a panel model for the first time to do a flow field calculation only (IRUN=0, EFLOW=l) 
and inspect the quality of the panel solution. A summary of the flow field calculation is contained 
on unit (FLSUM). Criteria such as smoothness of the pressure and vortex distributions are used to 
measure the quality of the flow field solution. Because of the relatively small execution time 
required for the flow field generation (90 seconds on the Cray XMP) and the ensured quality of the 
panel model, the flow field was generated in a single run along with the trajectory and ice accretion 
calculation (IRUN=5, IFLOW=l) and was stored (unit FLOWF). The job control input file (unit 
INPUT), the DUGLIFT geometry file (unit OGEOM), and the job summary file (unit JOBSUM) 
are shown in figures 24, 25, and 26 respectively. Figure 27 shows the pressure distribution at the 
0% spanwise location 

The second through fourth steps involve various steps in finding the ice accretion at each 
station of interest. These steps are repeated for each station of interest. 

The second step involves calculation of the streamline at the current section of interest. This 
requires finding the local stagnation zone, integrating upper and lower off-body streamlines from 
this stagnation zone and finding the on-body projection of the off-body streamline. This calculation 
required about 10 seconds for each section of interest. The coordinates for the off-body and on - 
body streamlines along with other information are contained in the job summary output (unit JOB- 
SUM). Figure 28 illustrates the off-body and on-body streamlines at the 0% span location. 

The collection efficiency at the station of interest is determined in the third step. This 
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involves determining the upstream release points of the droplets that would pass through the points 
at the leading edge demarcating the section of interest, determining the tangent trajectories associ- 
ated with these upstream release points, calculating the impact trajectories between these tangent 
trajectories, and calculating the collection efficiency resulting from these impact trajectories. The 
job s ummar y output file shown in figure 26 summarizes the pertinent information from these runs. 
Approximately 448 seconds on the Cray XMP were required to complete the trajectory calcula- 
tions for each section. Figure 29 depicts the impact trajectories generated in the collection effi- 
ciency calculation. 

During the fourth step the surface collection efficiencies generated in step three are used to 
find the collection efficiencies along the streamline generated in step two. Depending upon the 
value of NPSEC chosen, the determination of the collection efficiency values along the streamline 
are calculated from either an extrapolation (NPSEC=2) or an interpolation (NPSEC=4). In the cur- 
rent example very little spanwise variation was expected hence a “quasi-2D” collection efficiency 
calculation was made (NPSEC=2). The interpolation or extrapolation proceeds quickly and only 
required .1 seconds on the Cray XMP.The results from this calculation are summarized in the job 
summary output file (unit JOB SUM) in figure 26. Steps two through four are repeated for each of 
the sections of interest. 

Ice accretion along the streamline is calculated in step five for each of the streamlines. This 
involves calculating the ice thickness as a function of surface distance along the streamline at the 
current section of interest. The ice thickness distribution is calculated using a 3D version of the 
LEWICE heat transfer subroutine (ref. 7). The ice accretion at the local section is then calculated 
by adding the ice thickness calculated at a point to the point in either the surface normal or trajec- 
tory tangent direction.This results in a new off-body “iced-streamline” shown in figure 30. In addi- 
tion to the calculation of the “iced-streamline” the ice thickness distribution for each streamline is 
stored for the geometry modification calculation in step six. The ice accretion calculation, which 
proceeds fairly quickly (.2 seconds), is summarized in the job summary output file (unit JOBSUM) 
in figure 26. Step 5 is repeated for each of the sections of interest. 

The sixth and final step involves calculating the new goemetry from the ice thickness dis- 
tribution at each section of interest. This involves a cubic chordwise interpolation and a linear span- 
wise interpolation. The geometry modification calculations took approximately 1 second on the 
Cray XMP. Figure 31 shows the resulting iced wing resulting from the calculations. Figure 32 
shows a comparsion between the calculation and experiment for this case. In general the aggree- 
ment is good with the calculation predicting the shape, amount and postion of the ice. 
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Figure 22. - Swept NACA 0012 panel model used in example case. 
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COLLECTION EFFICIENCY 



SURFACE DISTANCE FROM HIGHLIGHT (CM) 


Figure 23. - Collection efficiency as a function of surface distance and spanwise location for 
the example case. 
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N012 

SWEPT NACA001 2 AIRFOIL (30 DEGREE SWEEP) EXAMPLE 
&TRAJ DFINE=.0001,DMDS(1)=20., 

VINF=67.0,CHORD=.4399, 

RHO=1 .29,RW=.1 ,TOL=.02,FNR=0.0, 

HMINI=.0001,HI=.01 > XSTART=-20.,XFINAL=1.5,TPRINT=.1,EPSI(1)=2.E-07 > 
EPSI(2)=2.E-07,EPSI(3)=2.E-07 &END 
&IMPING IRUN=5,ICE=1 ,ICEC=3,IFLOW=1 ,ISTRF=0,NBC(1 )=1 ,NBR(1 )=20, 
NBC(2)=1 ,NBR(2)=20,NBC(3)=1,NBR(3)=20,NPSEC(1)=2,NPSEC(2)=2, 
NPSEC(3)=2,SHFTF=0.,DSHIFT=.002,IRES=0, 

XSEC(1 ,1 )=0.,XSEC(1 ,2)=0.,XSEC{1 ,3)=0..XSEC(1 ,4)=0., 

YSEC(1 ,1 )=0.,YSEC(1 ,2)=0.,YSEC(1 ,3)=0.01 ,YSEC(1 ,4)=0.01 , 

ZSEC(1 ,1)=0.07,ZSEC(1 ,2)=-0.05,ZSEC(1 ,3)=0.07,ZSEC(1,4)=-0.05. 

XSEC(2,1)=0.318,XSEC(2,2)=0.318,XSEC(2,3)=0.,XSEC(2 > 4)=0. > 

YSEC(2,1)=0.55,YSEC(2,2)=0.55,YSEC(2 1 3)=0.01 1 YSEC(2,4)=0.01, 

ZSEC(2,1)=0.05,ZSEC(2,2)=-0.06,ZSEC(2,3)=0.07,ZSEC(2,4)=-0.05, 

XSEC(2,1)=-.318,XSEC(2,2)=-.318,XSEC(2,3)=0.,XSEC(2,4)=0., 

YSEC(2,1 )=-.55,YSEC(2,2)=-.55,YSEC(2,3)=0.01 ,YSEC(2,4)=0.01 , 
ZSEC(2,1)=0.05,ZSEC(2,2)=-0.06,ZSEC(2,3)=0.07 I ZSEC(2,4)=-0.05 1 
PLNST(1 ,1)=0.0,PLNST(1,2)=1 „PLNST(1 ,3)=0.,PLNST(1 ,4)=0., 
PLNST(2,1)=0.0,PLNST(2,2)=1.,PLNST(2,3)=0.,PLNST(2,4)=-.55, 

PLNST(3,1 )=0.0,PLNST(3,2)=1 .,PLNST(3.3)=0.,PLNST(3,4)=.55 &END 
&ICEIN LWC= .75 ,TAM B=267 . , PAM B=8987 6. , RH=1 00., 
XKINIT=.0045,SEGTOL=1.5,TSTOP=120.,DTIME=120.,DTFLOW=120. &END 


Figure 24. - Job control input file for example case unit (INPUT). 
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Figure 25. - Input flow field geometry for example case (unit OGEOM). 
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Figure 25. - Continued. 
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Figure 25. - Continued. 
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Figure 25. - Concluded. 
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Figure 26. • Continued. 
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Figure 26. - Continued. 



ooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooo 


o o o o o o o 


oooooooooooooooooooooo 


ooooooooooooo 
ooooooooooooo 
ooooooooooooo 
JE ooooooooooooo 
ry OOOOOOOOOOOOO 
|_ OOOOOOOOOOOOO 

in 

PQ OOOOOOOOOOOOO 


inNNOC^rte3N>TN«HC300Nin^NK>'3-v3r^OH«^00* 
inoifirHrHHfOLnotMflvofoocooocvcoinNON^HO^os 
^^^U^P^<>O^CvlfOU^sOCO(>0<— l(NJr-iOC0vT00CN'TLnC0OLnCT' 

hOCJCsj^o^oNar-Nou^^^cNjfMwoo'Ooinro^O'voLno'^^ 

>irJ)r -! ■ |nnnnnooOOOOOg'Q v C^(M>O v CO(CCOCOOP l 

JJJ^^^^.-^^^r^r^.-N.-MOOOOOOOOOOOO 


hio^nnnow^oioon 
'ON'ONHO'rO'OOsTvT «-HO 
Nfflv0 0'Nf0C'OOK10"0'0 
coin®i-«CM>JCMCM<Mr , '*oi--r- 
®r'-vDvoinM-fOCvj^O'r-'Oa'' 
C'CMMT'CMT'a'O'O'COCOCOCO 

OOOOOOOOOOOOO 


o 

O' 

<r 

va- 

VO 

O' 

CM 

o 

O 

CM 

CM 

r- 

•—< 

in 

r- 

o 

O 

c 

O' 

O' 

r-» 

c 

o 

to 

in 

® 

® 

c> 

o 


o 

r» 

NO 

vT 

CM 

o 

®vom 

to 

CM 

CM 

CM 

CM 

CM 

CM 

H 



H 


OOOOOOOO 


ooooooooooooooooooooo 


HWONOCOOlftHCO'tf ^ 
cor^O'OiocvirooM-®imno 
HvOOMOHNOCONNCNfO 
o3-cvj®0'0'®M'«~<r--'a- 
l0>a''0r«-0 , 'OCSjNTN00'>T'>T0^ 

OODOOHHHHHNNH 

ooooooooooooo 


K>®inr^vy^.--i'rw^W'3-®o®.H®vovyvoiMovr ;: soo©oop 

r^Sc'^r^ncoNoO'JcsiNNHv£oinotriH'or2--(0\TooD 
OC'C'O'O'OCOONNvO^^l'I’JMfONNr^^OOOOrvOOO 
ONC'C'^^OC'^r'O'O'CMJ'^O'IJ'^OC'C'OO'OCT'ff'ff'OOO 
0'C , 'C'0'0'0'0'C'^tM3'0'0'^ff'0'ff‘0'0'0'0'^0'^0'0'00 0 
OOOOOOOOOOOOOOOOOOOOOOOOOOi-Hr-li-H 


inMNT0"0^l0HHOHI^>T 

WOvO'OeOHvO'T'OCOH^O 

O'OHM'OHinoinovoHh- 

^ooooHHNNrort'ff'f 

00'0 S 0'0'0'P'®'0'0'0'^0> 

0sq\0'0^CT‘<7'O'O'0 v 0'0 s >0 < '0 s 

OOOOOOOOOOOOO 
I I 1 I I I I i I i I > ■ 


r^r^M^^inmi^C'Kii^©c-oc>rocMffl<reo©^<MvororN ; ©o© 

HCOh.^-CMiO'oniftNiflif"3'<?"3' N ®<MOinofvJONin>rooo 

inccrOHOsTocONincOHrovocooNC'OfflO'O'Oinooo 

SHOOHHNNroroc^cinininin^NONO'O'OvO'ONfOooo 

ooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooo 

I I I I I I I I I I I I I I I 1 * * I • * • 1 


s : 

oc. 

H- 

</) 

z 

> 


fsisoimncJr>*'Tf*«-so®«'3‘^ 

vOCMinfVj®©'— llAr'-rHCvICVIOv 
r^incMSJoincoo'ODr^'T^fO 
OOO'^flD'O'ffNOeO'OlOH 

fs. r«s. \o no vo 'O so *o in in m in 

OOOOOOOOOOOOO 
3 0 0 0 0 

fill 


^o^O'OvOHCOOvonNiriNNHvTOsjNr^chO'or^mvrooo 

inOrOLTN^HCWNCNIIT'OOmNN'TONHfflO'HH^OOO 

rOiriCOvO^^Nr^ONOfOCO'X(^Mr^^L^^NU\NOvHNHOOO 

sfNOOHMtf lf)W'ONNCOCOO'(>OOOHHHHNNvOOOO 

O 000000000000000 *^r-»I— «>-*r-»r -<0000 

ooooooooooooooooooooooooooooo 


jr fo^vOHHfOHinoror^oin 
fV vOin'OONO^O'T'TOrflOH 
f— LnHHCMT'HNOCT^r-vTP 
i/i csJC-ar-jor-tnoJcMn^r^roo 
z CJNNHHHHOOOO'ff'CO 

X HHHHHi-IHHrlHOOO 

ooooooooooooo 


^OOC^<>NC'Ov3rt'T^COOWCOOH(>0''OSOO^'TvfOJvO 

’Inroi^rOHHrHroLTOoN^o-NO'vrNOirir^rONNO^Nin^© 

rtoc'roNiHin^roN(\i'00'X«HUiOH'j , vocoo^oHHNM^ 

csjcsic\Jh^rO'v3-s3-'TLnLn'Ovor-.r-r-®®®a'C , 'a'0'0'Oooooo 

oooo*oo’ooooooooooooooooo^^^^-Hp-< 
i i i i i i i i i i » i i i i i * i * • 1 1 1 1 1 * * * 1 


rtMOCNj^Nr^C'O'WMrvN'OHcO'TC'TOCOMort'aoooo 

Nio^osrotor*sW'ONOiO(7''ri7'K)ON'omNHorooNoooo 

OHHHHOCO'OCHCDlAHCOCHNCHffl'O'JNHOOOOO 

inimnimmri'j-'TccroroKiNNtNiHHHOOoooooooo 

ooooooooooooooooooooooooooooo 

ooooooooooooooooooooooooooooo 


^U>NrHON^lOO'0'fO®HCO(>Oina)vON^O'^^^ S S“J5®? 

NTin'OfO^or^cNjsoooovroc'JC^^ofOOKi'OocsJ^oincO'^^. 

S-CMOr-crOunvOCOCrvO^HfvJKlfOsTvrva-rOtOCMCM^OOOOOO'O 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOO 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOO 
, I I I I I ( I I I I I I I 1 1 1 I i 1 t > 1 1 1 1 


%o>^ooorsi®^i^fOfocs»c^ro>Tcvf^®^K)incaroro<Nj^r^cvifo 

0'®®c>i^oin^cvir-wp-irsjror«-ror^‘a'srrsi'a'^<>r^ro^®M^ 

r4 (^rs.f0HH«Hr0WNONr0f0iHe0(vi'3 , K)0'HOC , wN^NU|\D 

HNl(^HinONK)NHL'>OvTCON>OP'fO<M7'H^'ONffl^^OH'J 

0JCviCMf0K)r0>T'3"imr>N0N0'0r-r»*p^®®®C'O0'CN0'C^c^ooo 

ooooo ooooooooooooooooooooo^.^^-. 


XC 


•« 


xc 


UJ 


xc 


o 








_> 


UJ 


in 

E 

1— 



Of 

LLt 

K> 

Q£ 

1— 

_J 

fN- 

O 

© 

CL 

CM 

u. 

M 

E 

• 



o 

o 

1- 


o 

va- 

=3 



in 

CL 


E 


1- 


UJ 


3 


ce 


O 


H 

it 


E 

© 


Q 

QC 



E 

E 

(/> 

UJ 

t-t 

UJ 

> 

Z 

H 

0£ 


H 


1- 


»“ 

3 

© 


3 

CL 



O 

u 

UJ 


QC 


z 


03 


l-H 


3 


1- 

E 

© 


3 

QC 



o 

I- 



QC 

© 

X( 


CQ 

X 

xc 


3 


xc 


© 



inNNO'OHHi^Noovr'O 

NCONf^'0(OOfOO"^ H '®' 0 

(MO'Hf^^rvHHIflHO'CDvO 

IflLftOC.O'O^ONCOOHPOin 

r^ovoojcivr-vo'3-iofocvir-to 

NNHHOOOOOOOOO 


vO'3‘'TCO'lTlCMn^f , 'J<— ‘fOvO 
coo'Or-iN\T'00'(M , oinr-ro 
sr^cocococDcocooO'O'crK) 
HHHHHHHH hj-i CM <NJ fO 
ooooooooooooo 
ooooooooooooo 

ooooooooooooo 

I I I I I I I t I I t I 1 


lomrotoooofotoovovinio 

(M^ininr'-rooN'OcooHO 

tsi>a“imn’C'ror-ifflsoiooin<-« 

cmcmcjcmcvjcvicvi.— 

inininininininu-iLniALninin 

ooooooooooooo 


^-ir^r^Ln^-.r-®'0'j-r^'Ovo\T 

rs.r>H\ovor*<7'roo"orsJCOo 

NO'r-ir-’or-oHinH<?'CON 

OOinv3'sT'TirifsJK)invO®0 

r'-oincMo , 'f--LnvTro<Mf— .oo 

LAirivT'TrorOfOrotOKirOKlKi 


o^cvro>Tinvoi^®^o^cv;rovrinvor-®C'OrHCsjrO'j-inNDr^® 

vo o *o *o 'C'^i'or^r s *M > *Nh* r* r>-r^-r*>®®®®®®®®® 


t-4 hnki^uI'Onco ^.ot-iCNJro 

q H HHH 


104 
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Figure 26. - Continued. 
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******* OUTPUT FROM SUBROUTINE CONFAC FOR ICEC= 
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- 20,000000 0.532183 - 0.007697 
- 20.000000 0.532193 - 0.009921 
- 20.000000 0.532202 - 0.011199 
- 20.000000 0.532211 - 0.012868 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Continued 
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SUBROUTINE ARYTRJ INPUT FOR SLICE: 3 

ROW ( NBRC) COLUMN (NBCC) XCNBRC.NBCC) Y(NBRC.NBCC) Z(NBRC,NBCC) 
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Figure 26. - Continued. 



SUBROUTINE BETAC OUTPUT FOR SLICE: 3 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Continued. 
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Figure 26. - Concluded. 
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Figure 27. - Pressure distribution at 50% span location for example case at 0 and 8 degrees 
angle-of-attack. 
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Figure 29. - Illustration of impact trajectories for the example case. 
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Figure 31. - Iced wing panel model for the example case. 
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Figure 32. - Comparison of predicted and measured ice shape at the 0 % span location for the 
example case. 


139 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No . 0704-0188 


Public reporting burden lor this collection of information is estimated to average 1 hour per response, including the time for reviewing instructor^,! searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any asp® 0 ! of this 
ooUection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate tor Information Operations and Reports, 121 5 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188). Washington, DC 20503. 


1. AGENCY USE ONLY ( Leave blank) 


2. REPORT DATE 

December 1993 


3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 


4. TITLE AND SUBTITLE 


Users Manual for the NASA Lewis Three-Dimensional Ice Accretion Code 
(LEWICE 3D) 


6. AUTHOR(S) 

Colin S. Bidwell and Mark G. Potapczuk 


5. FUNDING NUMBERS 


WU-505-68-10 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESSES) 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 


S. PERFORMING ORGANIZATION 
REPORT NUMBER 


E-7847 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESSES) 

National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 


10. SPONSORING/M ON ITORING 
AGENCY REPORT NUMBER 


NASA TM- 105974 


11. SUPPLEMENTARY NOTES 

Responsible person, Colin S. Bidwell, (216) 433-3947. 


12a. DISTRIBUTION/A VAl LABILITY STATEMENT 


12b. DISTRIBUTION CODE 


Unclassified - Unlimited 
Subject Categories 02, 03 


13. ABSTRACT (Maximum 200 words) 

A description of the methodology, the algorithms, and the input and output data along with an example case, for the NASA Lewis 3D 
ice accretion code (LEWICE3D) has been produced. The manual has been designed to help the user understand the capabilities, the 
methodologies and the use of the code. The LEWICE3D code is a conglomeration of several codes for the purpose of calculating ice 
shapes on three-dimensional external surfaces. A three-dimensional external flow panel code is incorporated which has the capability 
of calculating flow about arbitrary 3D lifting and nonlifting bodies with external flow. A 4th order Runge-Kutta integration scheme is 
used to calculate arbitrary streamlines. An Adams type predictor-corrector trajectory integration scheme has been included to calculate 
arbitrary trajectories. Schemes for calculating tangent trajectories, collection efficiencies and concentration factors for arbitrary regions 
of interest for single droplets or droplet distributions have been incorporated. A LEWICE 2D based heat transfer algorithm can be used 
to calculate ice accretions along surface streamlines. A geometry modification scheme is incorporated which calculates the new 
geometry based on the ice accretions generated at each section of interest. The three-dimensional ice accretion calculation is based on 
the LEWICE 2D calculation. Both codes calculate the flow, pressure distribution, and collection efficiency distribution along surface 
streamlines. For both codes the heat transfer calculation is divided into two regions, one above the stagnation point and one below the 
stagnation point, and solved for each region assuming a flat plate with pressure distribution. Water is assumed to follow the surface 
streamlines, hence starting at the stagnation zone any water that is not frozen out at a control volume is assumed to run back into the 
next control volume. After the amount of frozen water at each control volume has been calculated the geometry is modified by adding 
the ice at each control volume in the surface normal direction. 


14. SUBJECT TERMS 

Aircraft icing; Ice accretion prediction 

15. NUMBER OF PAGES 

143 

16. PRICE CODE 

A07 

17. SECURITY CLASSIFICATION 

IS. SECURITY CLASSIFICATION 

19. SECURITY CLASSIFICATION 

20. LIMITATION OF ABSTRACT 

OF REPORT 

OF THIS PAGE 

OF ABSTRACT 


Unclassified 

Unclassified 

Unclassified 



NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. Z39-18 
298-102 
















